i. ’EY KNOT IIBE-ARY 

V'lx ?0S’^‘ ' VJATE SCHOOL 

iJTEAJy, oi-iil'CliinA 95943-B002 



l)i.'T.EY 
AVAL P 
:OJ!JTEEi 



NAVAL POSTGRADUATE SCHOOL 
Monterey , California 




THESIS 



A STUDY OF SECOND AND THIRD ORDER 
MODELS FOR THE TRACKING 
SUBSYSTEM OF A RADAR GUIDED MISSILE 

by 

John W. Williams 

>s> Jt/ 

June 1988 

Thesis Advisor; H. A. Titus 

Approved for public release; distribution is unlimited 




I 



i 



UNCLASSIFIED 



ICURITY CLASSIFICA'^ION OF PAGE 



REPORT DOCUMENTATION PAGE 



3 . REPORT SECURITY CLASSIF. CATION 

UNCLASSIFIED 



lb RESTRICTIVE MARKINGS 



3. SECURITY CLASSIFICATION AUTHORITY 



,b DECLASSIFICATION /DOWNGRADING SCHEDULE 



3 DISTRIBUTION -AVAILABILITY OF REPORT 

Approved for public release; 
distribution is unlimited 



PERFORMING ORGANIZATION REPORT NUMBER(S) 



5 MONITORING ORGANIZATION REPORT NUMBER(S) 



■a. NAME OF PERFORMING ORGANIZATION 

Naval Postgraduate School 



6b OFFICE SYMBOL 
(If applicable) 

62 



7a NAME OF MONITORING ORGANIZATION 

Naval Postgraduate School 



Ic. ADDRESS {City, State, and ZIP Code) 

Vlonterey, California 93943-5000 



7b ADDRESS (C/ty, State, and ZIP Code) 

Monterey, California 93943-5000 



'a. NAME OF FUNDING / SPONSORING 
ORGANIZATION 



8b OFFICE SYMBOL 
(If applicable) 



9, PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER 



k. ADDRESS (Oty, State, and ZIP Code) 



10 SOURCE OF FUNDING NUMBERS 



PROGRAM 


PROJECT 


TASK 


WORK UNIT 


ELEMENT NO 


NO. 


NO 


ACCESSION NO 



\.TIUE (Include Security Classification) ^ STUDY OF SECOND AND THIRD ORDER MODELS FOR THE 
FRACKING SUBSYSTEM OF A RADAR GUIDED MISSILE 



2 PERSONAL AUTHOR(S) 

WILLIAMS, John W. 



13b. TIME COVERED 


14 DATE OF REPORT (Year, Month. Day) 


15 PAGE COUNT 


FROM TO 


1988 June 


139 



3a. TYPE OF REPORT 

Master's Thesis 



6 SUPPLEMENTARY NOTATioNThe views expressed in this thesis are those of the author 
ind do not reflect the official policy or position of the Department of 
^fens^ of th® U. S , Government . 



7 


COSATI CODES 


18 SUBJECT TERMS {Continue on reverse if necessary and identify by block number) 


FIELD 


GROUP 


SUB-GROUP 


radar guided missile, tracking, missile modeling. 








Coriolis equation, Kalman filter, aircraft 








modeling 



9 ABSTRACT {Continue on reverse if necessary and identify by block number) 

Phis thesis is a study of missile and target parameters used in second and 
third order modeling of the tracking subsystem used in radar guided missiles. 
Guidance methods are analyzed to determine which method is optimum in a 
search for an ’’ideal” missile. Target parameters which have an effect on the 
nissile tracking system are analyzed and a target acceleration probability 
nodel is discussed. A two dimensional third order tracking model is simu- 
lated utilizing a Kalman filter for target parameter estimation and 
prediction. Linear second and third order tracking models are simulated and 
compared with the third order Kalman filter tracker. This thesis concludes 
that a proportional navigation guidance method, with a non linear third order 
tracking Kalman filter, is the better model. Benefits of using a non linear 
third order Kalman filter may not override the cost and complexity of 
Implementation of the model. 



20 DISTRIBUTION /AVAILABILITY OF ABSTRACT 

ra UNCLASSIFlED/l4i;jL;MlTED □ PTIC USERS 




21. ABSTRACT SECURITY CLASSIFICATION 

UNCLASSIFIED 



22b TELEPHONE (Include Area Code) 

408-646-2560 



iZc OFFICE SYMBOL 

62Ts 



DD fORM 1473, 84mar 



83 APR edition may be used until exhausted 
All other editions are obsolete 

i 



SECURITY CLASSIFICATION OF THIS PAGE 

O U.S. Government Printing Office: 1986 — $06*24. 

UNCLASSIFIED 



Approved for public release; distribution is unlimited 

A STUDY OF SECOND AND THIRD ORDER MODELS FOR THE 
TRACKING SUBSYSTEM OF A RADAR GUIDED MISSILE 



by 



John Walter /^illiams 
Captain, United States Marine Corps 
B.S., United States Naval Academy, 1978 



Submitted in partial fulfillment of the 
requirements for the degree of 

MASTER OF SCIENCE IN ELECTRICAL ENGINEERING 



from the 



NAVAL POSTGRADUATE SCHOOL 
June 1988 



ABSTRACT 



This thesis is a study of missile and target parameters 
used in second and third order modeling of the tracking 
subsystem used in radar guided missiles. Guidance methods 
are analyzed to determine which method is optimum in a 
search for an "ideal” missile. Target parameters which have 
an effect on the missile tracking system are analyzed and a 
target acceleration probability model is discussed. A two 
dimensional third order tracking model is simulated 
utilizing a Kalman filter for target parameter estimation 
and prediction. Linear second and third order tracking 
models are simulated and compared with the third order 
Kalman filter tracker. This thesis concludes that a 
proportional navigation guidance method, with a non linear 
third order tracking Kalman filter, is the better model. 
Benefits of using a non linear third order Kalman filter may 
not overide the cost and complexity of implementation of the 
model . 



iii 



TABLE OF CONTENTS 



I. INTRODUCTION 1 

II. THE IDEAL MISSILE 3 

A. HEAD-ON ASPECT 6 

1. OG Target Acceleration 43 

2. 3G Target Acceleration 44 

3. 6G Target Acceleration 45 

B. TAIL ASPECT 45 

1. OG Target Acceleration 45 

2. 3G Target Acceleration 46 

3. 6G Target Acceleration 46 

C. BEAM ASPECT 47 

1. OG Target Acceleration 47 

2. 3G Target Acceleration 48 

3. 6G Target Acceleration 48 

D. CONCLUSIONS 4 8 

III. TARGET MODEL 50 

A. TACTICS FOR MISSILE DEFENSE 53 

IV. SYSTEM MODEL 56 

A. COORDINATE SYSTEM 56 

B. EQUATIONS OF MOTION 57 

C. SECOND ORDER MODEL 60 

D. THIRD ORDER MODEL 62 

V. KALMAN FILTER 64 

VI. SIMULATION 70 

A. ASSUMPTIONS 70 

B. INITIALIZATION 71 

C. SECOND ORDER MODEL 71 

D. THIRD ORDER MODEL 71 

1. Time Invariant, Constant Phi Model . . 72 

2. Time Variant, Variable Phi Model ... 72 

E. RESULTS 7 5 

1. Gains 75 



IV 



2. Head-on Aspect 107 

a. OG Target Acceleration 107 

b. 3G Target Acceleration 107 

c. 6G Target Acceleration 108 

3. Tail Aspect 108 

a. OG Target Acceleration 108 

b. 3G Target Acceleration 108 

c. 6G Target Acceleration 109 

4. Beam Aspect 109 

a. OG Target Acceleration 109 

b. 3G Target Acceleration 109 

c. 6G Target Acceleration 110 

VII. CONCLUSIONS Ill 

REFERENCES 112 

APPENDIX A. IDEAL MISSILE PROGRAM LISTING 113 

APPENDIX B. THIRD ORDER SIMULATION PROGRAM 

LISTING 117 

INITIAL DISTRIBUTION LIST 130 



V 



TABLE OF VARIABLES 



AT 

AM 

B 

BO 

BDO 

BDOT 

BDDOT 

D2R 

DELR 

DXT 

DYT 

G 

GAMMA 

GO 

GDO 

GDOT 

GDDOT 

GR 

GS 

H 

K2F 

LOS 

LOSD 

LPHDG 

LPHDGO 

LPLOS 

LPLOSD 

LPR 

LPRD 

LPXM 

LPYM 

MISSXO 



acceleration of target (g) 
acceleration of missile (g) 
antenna angle 

initial antenna angle (rad) 

initial antenna angular rate (rad/sec) 

antenna angular rate 

antenna commanded angular acceleration 
conversion factor degrees to radians 
matrix of state space range errors 

X coordinate, impact point for Direct Path missile 

Y coordinate, impact point for Direct Path missile 

force of gravity, 32.2 ft/sec/sec 

antenna angle (rad) 

initial antenna angle (rad) 

initial antenna angular rate (rad/sec) 

antenna angle rate (rad/sec) 

antenna commanded angular acceleration 

matrix of gains, range channel for Kalman filter 

matrix of gains, bearing channel for Kalman filter 

observation vector, used in comment lines 

conversion factor knots to feet/sec 

line of sight angle (rad) 

line of sight rate (rad/sec) 

lead pursuit missile heading (rad) 

initial lead pursuit missile heading (rad) 

line of sight for lead pursuit missile (rad) 

line of sight rate for lead pursuit missile 

range for lead pursuit missile to target (ft) 

range rate for lead pursuit missile (ft/sec) 

lead pursuit missile position X coordinate 

lead pursuit missile position Y coordinate 

initial missile location X parameter (ft) 



VI 



MISSYO 


initial missile location Y parameter (ft) 


P 


prediction covariance matrix. 


PHDG 


initial proportional navigation missile heading 


PHI 


discrete model state space matrix 


PNAM 


prop nav missile commanded acceleration 


PNHDG 


prop nav missile heading (rad) 


PNLOS 


line of sight angle for prop nav missile (rad) 


PNLOSD 


line of sight rate for prop nav missile (rad/sec) 


PNR 


range for prop nav missile to target (ft) 


PNRD 


range rate for prop nav missile (ft/sec) 


PPHDG 


pure pursuit missile heading (rad) 


PPLOS 


line of sight angle for pure pursuit missile (rad) 


PPLOSD 


line of sight rate for pure pursuit missile 


PPR 


range for pure pursuit missile to target (ft) 


PPRD 


range rate for pure pursuit missile (ft/sec) 


PPXM 


pure pursuit missile position X coordinate 


PPYM 


pure pursuit missile position Y coordinate 


PR 


error covariance matrix for range channel 


PROP NAV 


proportional navigation 


PS 


error covariance matrix for bearing channel 


PXM 


prop nav missile position X coordinate 


PYM 


prop nav missile position Y coordinate 


QR 


maneuvering covariance matrix range channel 


QS 


maneuvering covariance matrix bearing channel 


RO 


initial range missile to target (ft) 


R,RM 


measured range from missile to target (ft) 


RD 


range rate from missile to target (ft/sec) 


RDDK 


estimate of range acceleration at time k 


RDDKPl 


prediction of range acceleration for time k+1 


RDK 


estimate of range rate at time k 


RDOTM 


measured range rate 


RKPl 


prediction of range rate for time k+1 


RK 


estimate of range at time k 


RKPl 


prediction of range at time k+1 


RMCOV 


measurement covariance matrix for range channel 



Vll 



RNG 

RPHI 

S 

S 

SDEL 

SG 

SMCOV 

SOHDG 

SOLOS 

SOR 

SOU 

SOXM 

SOYM 

TEMPI 

TEMP2 

TEMP3 

TGTA 

TGTAX 

TGTAY 

TGTHDG 

TGTV 

TGTVXO 

TGTVYO 

THDG 

TK 

TK 

TTGO 

TTG 

TVELX 

TVELY 

U 

VM 

VT 

XT 

YT 



matrix of range state space values 
PHI matrix for range channel 
Sigma, estimate of line of sight 
matrix of state space sigma values 
matrix of error values 

gain matrix for bearing channel for Kalman filter 
measurement covariance matrix, bearing channel 
missile heading of second order simulation (rad) 
line of sight for second order simulation (rad) 
range for second order simulation (ft) 
commanded acceleration for second order simulation 
missile position X coordinate for second order 

missile position Y coordinate for second order 
temporary matrix for storing values 
temporary matrix for storing values 
temporary matrix for storing values 
acceleration of target (ft/sec/sec) 
target acceleration in X direction 
target acceleration in Y direction 
target heading (rad) 
velocity of target, 667 ft/sec 

initial velocity of target X parameter (ft/sec) 

initial velocity of target Y parameter (ft/sec) 

initial target heading (rad) 

sample interval (sec) 

square of the sample interval 

initial time to go (sec) 

time to go until impact (sec) 

target velocity in X direction (ft/sec) 

target velocity in Y direction (ft/sec) 

commanded missile acceleration 

velocity of missile, 2500 ft/sec 

velocity of target (kts) 

target position X coordinate 

target position Y coordinate 



Vlll 



INTRODUCTION 



I . 



Aviation plays an extensive role in current combat 
scenarios. An aircraft, because of its capability to carry 
missiles, is a very formidable weapon platform. Missiles 
provide offensive killing power which change tactics in 
battle scenarios. In order to have the edge in the air to 
air arena, an aircraft must possess the best defensive and 
offensive capabilities. One of the main weapons in the 
aircraft arsenal is the radar guided missile. 

Radar guided missiles are all weather capable, can be 
launched outside of visual range and are less susceptible to 
countermeasures, compared with other missile types. 
Improvement of the missile is a constant necessity to 
maintain air superiority. 

Improvements in aircraft maneuverability dictate the 
need for missiles to increase performance and capabilities. 
A rule of thumb for design is that the missile must have a 
4:1 acceleration advantage over the target. With modern 
aircraft able to sustain lateral accelerations of ten times 
the force of gravity (G) the missile must be capable of 40G. 

Guidance methods are chosen which optimize the missile 
capabilities to destroy the target. The guidance method 
selected has a large impact on the design of other missile 
subsystems. For any missile to guide to the target, the 
sensor subsystem must track the target. To optimize the 
guidance and tracking of radar guided missiles a predicting 
filter can be used. 

One of the simplest missile guidance techniques is a 
second order system which compensates for bearing error and 
bearing rate error. This thesis will look at third order 
models to help optimize the missile sensor subsystem to 
provide better guidance command inputs. A major impetus for 
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finding an optimum predicting guidance method is to improve 
missile performance in the final guidance stages. 

There is a region at the end of the missile flight path, 
where the time to intercept is so short that inputs to the 
missile control surfaces will not be effective. If the 
missile has an exact solution of target parameters, it can 
predict the future target position, through the time where 
no inputs will be effective. The missile is flying to the 
projected position of the target, compensating for the time 
delay of control effectiveness. 

Section II will look at some guidance methods and target 
parameters to aid in finding the optimum missile guidance. 
Section III looks at target models and how to implement them 
in the computer simulation. Section IV derives the second 
and third order models. A Kalman Predicting Filter is 
discussed in Section V. Computer simulation and 
implementation of two dimensional third order models with 
the Kalman Filter is given in Section VI. Section VII 
contains conclusions and recommendations. Program listings 
of the computer simulation are in Appendix A and Appendix B. 
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II. THE IDEAL MISSILE 



Definition of an ideal missile is very difficult. Cost 
or performance functions can be generated to account for 
miss distance, fuel expended, control inputs, flight path 
and numerous other parameters. To the operator the ideal 
missile is one that destroys the target. The designer must 
try to account for a variety of scenarios and targets to 
design the optimum missile. Tradeoffs of performance and 
costs will dictate what the final sub-optimum missile will 
be . 

In an attempt to find an ideal missile insight may be 
gained by looking at the various flight paths and parameters 
for different guidance methods. Three methods to look at 
are pure pursuit, lead pursuit and proportional navigation. 
Pure pursuit would entail a missile always flying directly 
at the target. Lead pursuit would be the case of a missile 
always flying to a point slightly ahead of the target. The 
magnitude of the lead may vary from one time constant, to 
the total missile flight time, ahead of the target. The 
latter would produce a lead collision intercept. 
Proportional navigation utilizes line of sight rate to guide 
the turn of the missile, to zero any further line of sight 
rate. The commanded acceleration to the missile is given by 
the equation: 

AM = k p Vc (1) 

where AM = missile acceleration 

k = constant of proportionality 
P = antenna angle rate 
Vc = closing velocity 

The constant of proportionality is determined by the 
designer. Observing the effects of the proportionality 
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constant on the line of sight rate shown in Figure II-l, any 
constant above 3 is an appropriate value. For k less than 3 
the line of sight rate has its largest slope at the end of 
the intercept. For k greater than 3 the line of sight rate 
will be very small prior to the impact point. 




Figure II-l. Prop Nav Proportionality Constant [Ref. 1] 



Three guidance methods are analyzed to compare flight 
paths, heading changes, line of sight angle and line of 
sight rate. A fourth missile is simulated and plotted, 
identified by "direct". 

The direct path missile is programmed using "a 
posteriori" knowledge of the target flight path. The direct 
path missile goes to the point of impact, in a straight 
line, from the point of launch. 

If the "ideal missile" is defined such that there is no 
missile maneuvering, burns minimal fuel, has the greatest 
launch distance, minimum intercept time and accounts for all 
target maneuver, it will be a direct path missile. 
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Two programs were written using Dynamic Simulation 
Language, DSL, to analyze the guidance methods and produce 
plots. The programs are included in Appendix A. 

Three scenarios are used to compare the guidance 
methods : 



head-on aspect 
tail aspect 
beam aspect 

Three target accelerations are used to determine the 
effect of the target maneuver on the parameters. The 

selected target accelerations are: 

0 g ' s 
3 g ' s 
- 6 g's 

A second order proportional navigation missile with a 
proportional navigation guidance constant of four, where 
the line of sight angle and angle rate is estimated by the 
antenna parameters of angular position and angular rate. 
The governing equations for tracking in bearing are: 



P 



'T . 

p dt + po 
o 



( 2 ) 



P 



‘T . . 

P dt + 3o 
. o 



(3) 



P'= -20*B + 100*(LOS-P) 



(4) 



where 



P = antenna position angle 
P = antenna angle rate 

P = antenna angular acceleration 

LOS = actual target bearing 



Pure pursuit and lead pursuit guidance missiles are 
initialized heading directly at the target. Pure pursuit 
guidance maintains heading directly at the target by 
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calculating 


the heading at each step of the discrete 


simulation . 


The equation for pure pursuit heading is: 


PPHDG = 


atan ( (yt-ym) / (xt-xm) ) (5) 


where 


yt = current target Y coordinate 
xt = current target X coordinate 
ym = current missile Y coordinate 
xm = current missile X coordinate 



Lead pursuit maintains a heading in front of the target, 
with a variable lead, calculated using half the time to go. 



given by: 




LPHDG = 


atan( (yt+tvely*.5ttg)/(xt+tvelx*.5*ttg) ) (6) 


where 


yt = current target Y coordinate 

xt = current target X coordinate 

tvely = target velocity in Y direction 
tvelx = target velocity in X direction 
ttg = time to go in the intercept 



The results of the comparisons are given in graph form 
in Figure II-2 through Figure 11-37. 



A. HEAD-ON 


ASPECT 


Figures 


II-2 through 11-13 are the results of missile 



guidance comparisons for head-on aspect initial condition 
with OG, 3G and 6G constant target acceleration. The 
missile begins at the origin of the graph. The target 
initial position is at x = 10000 ft, y = 500 ft. Applied 
lateral target acceleration is away from the missile. 
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Figure II-3 



Missile Heading Plot for Head-on 
Aspect No Target Turn 
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Figure II-4 Line of Sight Angle Plot for Head-on 
Aspect No Target Turn 
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Figure II-5 Line of Sight Rate Plot for Head-on 
Aspect No Target Turn 



C 







(iJ) Jl 



Figure 11-6 
Aspect 



Position Plot for Head-on 
Constant 3G Target Turn 
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Figure 



II-7 Missile Headina Plot for Head-on 
Aspect Constant 3G Target Turn 
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Figure II-8 Line of Sight Angle Plot for Head-on 
Aspect Constant 3G Target Turn 
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Figure II-9 Line of Sight Rate Plot for Head-on 
Aspect Constant 3G Target Turn 
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Figure II-IO Position Plot for Head-on 
Aspect Constant 6G Target Turn 
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Figure II-ll Missile Headings Plot for Head-on 
Aspect Constant 6G Target Turn 
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Figure 11-12 Line of Sight Angle Plot for Head-on 
Aspect Constant 6G Target Turn 
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Figure 11-13 Line of Sight Rate Plot for Head- 
Aspect Constant 6G Target Turn 
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Figure 11-14 Position Plot tor Tail 
Aspect Constant No Target Turn 
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Figure 11-16 Line of Sight Angle Plot Tail 
Aspect Constant No Target Turn 
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Figure 11-17 Line of 
Aspect Constant 



Sight Rate Plot Tail 
No Target Turn 
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Figure 11-18 Position Plot Tail 
Aspect Constant 3G Target Turn 
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Figure 11-19 Missile Heading 
Aspect Constant 3G Target 



Plot Tail 
Turn 
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Figure 11-20 Line of 
Aspect Constant 



Sight Angle Plot Tail 
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Figure 11-21 Line of Siaht Rate Plot Tail 
Aspect Constant 3G’Target Turn 
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Figure 11-22 Position Plot Tail 
Aspect Constant 6G Target Turn 
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Figure 11-23 Missile Heading Plot Tail 
Aspect Constant 6G Target Turn 
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Figure 11-24 
Aspect 



Line of 
Constant 



Sight Angle Plot Tail 
6G Target Turn 
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Figure 11-25 Line of Siaht Rate Plot Tail 
Aspect Constant EG’Target Turn 
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Figure 11-26 Position Plot Beain 
Aspect Constant No Target Turn 
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Figure 11-28 Line of Sight Angle Plot Beam 
Aspect Constant No Target Turn 
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Figure 11-29 Line of Sight Rate Plot Beam 
Aspect Constant No Target Turn 



34 



o O 




(ij) ; 



Figure 11-30 Position Plot Beam 
Aspect Constant 3G Target Turn 
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Figure 11-31 Missile Heading 
Aspect Constant 3G Target 
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Figure 11-32 Line of Sight Angle Plot Be 
Aspect Constant 3G Target Turn 
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Figure 11-33 
Aspect 



Line of Sight Rate Plot Beam 
Constant 3G Target Turn 
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Figure 11-34 Position Plot Beam 
Aspect Constant 6G Target Turn 
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Figure 11-35 Missile Headina Plot Beam 
Aspect Constant 6G Target Turn 
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Figure IT-36 Line of Sight Angle Plot Beair. 
Aspect Constant 6G Target Turn 
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Figure 11-37 Line of Sight Rate Plot Beam 
Aspect Constant 6G Target Turn 
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By applying the acceleration away from the missile the 
pilot will lose sight of the missile. This is undesirable, 
but a turn into the missile will help the missile, in the 
early stages, more than a turn away. With a turn into the 
missile, the pilot will also lose sight of the missile, 
during a constant acceleration turn. 

1 . OG Target Acceleration 

From the position plot. Figure II-2, it is seen that 
proportional navigation guidance is essentially the same as 
the direct path guidance. Any errors are due to 
initialization of the heading for the proportional 
navigation guidance. The pure pursuit guidance missile and 
the lead pursuit guidance missile fly curvilinear paths to 
targetintercept . The curve for the lead pursuit guidance 
missile is less than the pure pursuit missile due to target 
lead . 

The missile headings graph. Figure II-3, shows the 
relative heading changes involved for each missile guidance. 
After the initialization errors have been corrected, the 
proportional navigation guidance missile parallels the 
direct path missile. The heading changes for the lead 
pursuit are less than for the pure pursuit guidance method. 
Large increases in missile headings at the end of the 
intercept implies large lateral accelerations are required 
for the missile to complete the intercept. 

The line of sight graph. Figure II-4, shows what 
would be expected for this case. The proportional 
navigation guidance and direct path missile maintain 
constant line of sight, approximately, while the line of 
sight increases for lead pursuit and pure pursuit guidance 
methods. The large change in line of sight at the end of 
the intercept also correlates to a high lateral acceleration 
required by the missile. 

The line of sight rate graph. Figure II-5, gives 
some insight to the control inputs to the missile guidance 



43 



M-rVvoouCtL/ Ml OUVtHNMENT EXPENSE 



subsystem. The values from, Figure II-5, are the slew rated 
for the sensor subsystem. A positive slew rate is seen for 
the prop nav guidance only at the final stages of the 
intercept. Lead pursuit and pure pursuit guidance methods 
have accelerating positive slew rates throughout the 
intercept. The direct path missile has a negative slew 
rate, caused by the missile speed advantage (2500 : 667 
ft/sec) . 

2 . 3G Target Acceleration 

The position graph, Figure II-6, shows a curvilinear 
path for all three guidance methods. The curvature of the 
target flight path is misleading, because of the axis 
scaling. The target is maintaining a constant acceleration. 
All three guidance methods appear to end up in a tail chase. 
The scaling is misleading again. Target heading change is 
approximately 50 degrees. The proportional navigation 
missile impacts in the beam while lead pursuit will be rear 
quarter and pure pursuit will be a tail chase. 

The missile headings graph. Figure II-7, shows the 
proportional navigation guidance has the lowest heading 
slope and is approximately linear at the end of the 
intercept. Pure pursuit and lead pursuit guidance methods 
have accelerating missile heading slopes requiring higher 
missile acceleration. 

The line of sight graph. Figure II-8, is similar to 
the Figure II-4, proportional navigation guidance method, 
which has low line of sight angles, slightly increasing due 
to target acceleration. Lead pursuit and pure pursuit 
guidance methods have line of sight angles which increase at 
an accelerated rate throughout the intercept. The direct 
missile has decreasing line of sight. The large negative 
LOS for the direct missile at the end of the intercepts is 
caused by the miss distance and heading initialization. 

Line of sight rates. Figure II-9, correlate with the 
line of sight plot. Figure II-8. Line of sight rates are 
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small but show an acceleration at the end of the intercept, 
due to decreasing range. Proportional navigation guidance 
methods are reducing the line of sight angle while lead 
pursuit and pure pursuit increase the line of sight angle. 

3 . 6G Target Acceleration 

Comparing the position plot, Figure II-IO, with that 
of the 3G case. Figure II-6, similar statements can be made 
about all of the missile paths. Scaling is slightly 
deceiving; the target has made approximately 100® heading 
change. All flight paths are curvilinear with the 
proportional navigation guidance method being the shorter of 
the three methods. 

Figure II-ll, shows the heading changes for the 
missiles and the smaller missile maneuvering required for 
the proportional navigation missile. Figure 11-12 and 
Figure 11-13 show larger magnitudes for line of sight angle 
and line of sight rate than the 3G case, but follow the same 
trends. The direct path missile shows a reversal in line of 
sight rate as the target heading change is greater than 90°. 

B. TAIL ASPECT 

Figure 11-14 through Figure 11-25 are the results of 
missile guidance comparisons for tail aspect initial 
condition with OG, 3G and 6G constant target acceleration. 
The missile begins at the origin of the graph. The target 
initial position is X=10000 ft. and Y=1000 ft., with an 
initial heading of 090, parallel to the X axis. Applied 
target acceleration is directed into the missile, 
perpendicular to the target heading. 

1 . OG Target Acceleration 

The position plot, Figure 11-14, shows the 
proportional navigation guidance missile flies a similar 
path as the direct path missile. The difference in the 
flight paths is due to errors in initialization. The 
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missile heading plot. Figure 11-15, shows that the heading 
for proportional navigation guidance and direct path 
missiles are parallel after the initialization errors are 
corrected. Pure pursuit and lead pursuit guidance methods 
have continually changing headings with accelerating slopes 
at the final stage of the intercept. 

Line of sight angles for the proportional navigation 
guidance and direct path are approximately constant and 
equal to the initial line of sight angle, giving a constant 
bearing decreasing range trajectory as seen by the target. 
The line of sight angle for pure pursuit and lead pursuit 
guidance decrease, but non linearly, as seen in Figure 11-16 
and Figure 11-17 . 

2 . 3G Target Acceleration 

With target acceleration, all three missiles fly a 
curvilinear path. The proportional navigation guidance 
method has the shortest flight path, as seen in Figure Il- 
ls. Proportional navigation guidance gives a linear heading 
change, as seen in Figure 11-19. Pure pursuit and lead 
pursuit guidance methods give higher heading slopes when the 
target applies lateral acceleration as compared to the OG 
heading plot. Figure 11-15. 

Line of sight angle changes are small for propor- 
tional navigation guidance, as shown in Figure 11-20. Pure 
pursuit and lead pursuit guidance have decreasing line of 
sight angles, with corresponding decreasing line of sight 
rates, as seen in Figure 11-20 and Figure 11-21. Line of 
sight rates increase for proportional navigation guidance, 
as would be expected from the path the missile flies. The 
direct path gives both a nonlinear line of sight angle and 
line of sight rate throughout the intercept. 

3 . 6G Target Acceleration 

When the target acceleration is increased, flight 
paths have a larger curvature, as seen in Figure 11-22. The 
scaling gives some distortion, the target has gone through 
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approximately 100° of heading change. The missile heading 
changes, as per Figure 11-23, are similar to those observed 
for the 6G head-on aspect. Figure II-ll. 

The line of sight angle and line of sight rate are 
larger for an increase in lateral target acceleration, as 
seen by comparing Figure 11-24 and Figure 11-25 with the 3G 
case, Figure 11-20 and Figure 11-21. The line of sight and 
line of sight rate for the direct path have a very large 
slope at the final intercept due to effects of decreased 
range. The reversal of line of sight rate for the direct 
path in Figure 11-25 is where the target heading change is 
90° . 

C. BEAM ASPECT 

Figures 11-26 through 11-37 are the result of missile 
guidance comparisons for beam aspect initial conditions with 
OG, 3G and 6G constant target accelerations. The missile 
begins at the origin of the graph. The target initial 
position is X=15000, Y=0 . Applied acceleration is directed 
into the missile. 

1 . O G Target Acceleration 

As in the two previous cases, with no lateral target 
acceleration, proportional navigation guidance and direct 
path missiles have similar flight paths. Pure pursuit and 
lead pursuit guidance have curvilinear flight paths, as seen 
in Figure 11-26. Heading changes are small for proportional 
navigation guidance and direct flight path missiles but not 
zero as what might be inferred from Figure 11-27, because of 
scaling. The heading change for pure pursuit and lead 
pursuit guidance is accelerating throughout the flight time 
with the intercept ending in a tail chase. 

Line of sight and line of sight rate. Figure 11-28 
and Figure 11-29, are similar to the two previous cases, for 
no target acceleration and the analysis is the same. 
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2 . 3G Target Acceleration 

For the beam aspect initial condition, when lateral 
acceleration is applied, the difference between flight paths 
for proportional navigation guidance and direct path is 
opposite from the two previous cases for lateral 
acceleration. The proportional navigation guidance missile 
flight path is on the opposite side of the direct path from 
pure pursuit and lead pursuit guidance flight paths. Figure 
11-30, with the opposite curvature. Headings for 
proportional navigation guidance continually decrease while 
pure pursuit and lead pursuit guidance increase. 

Differences between the methods are enhanced by the 
line of sight angle and line of sight rate plots in Figure 
11-32 and Figure 11-33. Proportional navigation guidance 
decreases line of sight while the others have an increasing 
line of sight and appropriate line of sight rate. 

3 . 6G Target Acceleration 

Increasing the target lateral acceleration magnifies 
the flight path differences between the guidance methods. 
As has been seen from the previous cases, the larger lateral 
acceleration increases the magnitudes of the values for 
Figure 11-33 through Figure 11-37, compared with similar 
graphs from the other cases, but the trends remain the same. 
Proportional navigation guidance parameters have smaller 
changes than pure pursuit and lead pursuit guidance, with 
parameters generally decreasing instead of increasing for 
pure pursuit and lead pursuit guidance. 

D. CONCLUSIONS 

For scenarios with no applied target lateral 
acceleration the proportional navigation missile is the same 
as the direct missile. The pure pursuit and lead pursuit 
missiles finish in a 'tail chase’ where a missile speed 
advantage is required to complete the intercept. The line 
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of sight remains constant for the proportional navigation 
missile but increases with the pursuit missiles. The line 
of sight rate increases with decreased range for the pursuit 
missiles but is zero for the proportional navigation and 
direct missiles. 

When target lateral acceleration is applied, there is a 
deviation between the direct and the proportional navigation 
missiles. Since the target is turning into the missile, the 
line of sight angle decreases at an accelerated rate as 
range decreases. The proportional navigation missile 
accounts for the change of line of sight by turning into the 
target. The pursuit missiles fly a tail chase profile with 
higher line of sight accelerations due to the target turn. 

For the direct missile, when target acceleration is 
applied, the line of sight is not constant and the rate of 
change depends on the applied acceleration. Implementation 
of a direct missile is impossible because the parameters 
used to guide the missile are dependent on the target flight 
path . 

If an optimum missile is to be designed, proportional 
navigation guidance is the closest to an "ideal" missile. 
The better the proportional navigation missile can 
compensate for the effects of the target acceleration, the 
closer to "ideal" the missile will become. 
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Ill , TARGET MODEL 



A complete target model for use in computer simulation 
is very involved, time consuming and computer intensive. To 
simplify target simulation the target flight profile is 
based on the fact that the missile sees only the effects of 
the target command inputs and resultant flight path. The 
target model was simplified to include only the flight 
profile desired and not be concerned with the full target 
modeling. A constant speed, constant acceleration target is 
assumed for the simulation. A variable speed, variable 
acceleration target can be added at a later time. The 
missile simulation estimates and predicts target parameters 
of range, range rate, range acceleration, bearing, bearing 
rate, and bearing acceleration. Therefore, for proper 
evaluation of the missile guidance and missile flight 
profiles the target parameters in missile coordinates for 
all these parameters must be computed. 

Complete analysis of target motion is obtained from a 
three dimensional derivation, but insight can be gained from 
two dimensional modeling. Three dimensional flight profiles 
are easily implemented on the computer but graphic display 
of the results are difficult. Two dimensional displays are 
easier to implement and comprehend. A two dimensional 
target model is assumed. 

A target can accelerate at values ranging from negative 
maximum instantaneous acceleration to positive maximum 
instantaneous acceleration, A(max inst) . A(max inst) is 
defined as the aerodynamic acceleration given by the maximum 
deflection of control surfaces. A (max inst) is dependent on 
airspeed and air density. High speeds and low altitudes 
produce the highest instantaneous accelerations. Maximum 
sustained acceleration, A(max sust) , is defined as the 
aerodynamic acceleration to maintain constant airspeed and 
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constant altitude, at full thrust. A(max sust) can be 
exceeded but must be compensated for by a reduction in 
airspeed or altitude. 

In three dimensional maneuvering cross coupling exists 
between horizontal and vertical angle and angle rate 
components of target velocity and target acceleration. 
Applied accelerations and velocity changes in one direction 
will affect the parameters, seen by a missile, in the two 
other directions. 

Thrust capabilities have a direct correlation to A(roax 
sust) and the airspeed of an aircraft. An aircraft with 
higher thrust can maintain a higher speed and compensate for 
drag induced by the applied acceleration. A modern aircraft 
with a relatively high thrust to weight ratio will have a 
very high A (max sust) which is close to A (max inst) . 

Airspeed is a key element for maneuverability and 
survivability. Tactics incorporate optimum techniques for 
maintaining airspeed or recovering lost airspeed. Pilots 
learn to compensate for limitations of A (max sust) by 
intentionally decreasing altitude and use the effects of 
gravity to maintain airspeed when lateral acceleration is 
applied. Another technique is to apply the lateral 
acceleration required to perform a maneuver then to reduce 
the acceleration, allowing excess thrust to restore the 
airspeed and altitude lost during the maneuver. There is a 
recovery time for the thrust to restore the lost energy, so 
to aid in restoring airspeed, a pilot will normally go to 
zero acceleration, reducing any induced drag, effectively 
increasing the aircraft thrust. This maneuver causes a loss 
in altitude due to gravity but improves airspeed 
restoration . 

The probability of aircraft acceleration is used to 
determine parameters for the target model used in missile 
designs and simulations. Figure III-l shows a typical 
acceleration probability graph used in missile design. The 
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figure does not account for pilot tendencies nor the 
difference between A(niax sust) and A{max inst) . The graph 
assumes that a pilot will maneuver primarily at zero 
acceleration, straight and level, or maximum acceleration, 
for a turn, with some probability for any other possible 
acceleration. The design engineer assigns probabilities for 
the impulse functions at zero acceleration and at maximum 
acceleration depending on the type of target aircraft. A 
large bomber may have an A (max sust) half that of a fighter 
aircraft, with less probability of turning than flying 
straight and level. 
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A. TACTICS FOR MISSILE DEFENSE 



When a missile is fired at the target the battlefield 
scenario changes to an immediate survival situation. If the 
missile is undetected, the acceleration probabilities of 
Figure III-l may be an adequate target model. If a pilot 
sees the missile, pilot reactions will change the 
probabilities. How the probabilities change may be of 
consequence to the missile guidance. An optimum missile 
design may be able to use pilot tendencies to increase 
missile performance. The overall acceleration probability 
from Figure III-l is a zero mean function with a variance, 
o2, dependent on the probabilities assigned. Reference 2 
discusses obtaining the parameters for the target accel- 
eration probability model of Figure III-l. 

When the pilot imposes a missile defense, the overall 
acceleration may not be zero mean, nor maintain the same 
variance. To account for the changes in acceleration 
probability a function similar to Figure III-2 might be 
used. This model accounts for some variance to the 
acceleration which the pilot is trying to achieve centered 
around A (max sust) . If the pilot is trying to achieve A (max 
inst) , it is assumed he will be decreasing airspeed and 
reducing actual acceleration until the applied acceleration 
is decreased to A(max sust) or below. A smart pilot will 
either fly at a maximum acceleration or at zero 
acceleration, increasing the aircraft maneuverability. 

Last ditch maneuvers are performed at A (max inst) to 
avoid the missile, neglecting any adverse effects of 
applying acceleration, in order to increase survivability. 
If the last ditch maneuver is performed too soon, 
acceleration is decreased, due to loss of airspeed, negating 
the effectiveness of the maneuver. Further studies can be 
made correlating the use of A(max inst) versus A(max sust) 
for missile defense tactics. 
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Figure III-2 Probability of Acceleration 

If pilot reaction is taken as Gaussian when trying to 
achieve a desired acceleration, the overall acceleration 
probability will be Gaussian with a non-zero mean and a 
variance dependent on the combination of the two Gaussian 
terms. The probability assigned to each term determines the 
mean and variance. 

Missile defense includes placing the missile on the 
beam, to utilize the largest acceleration vector, with 
maneuvers made out of phase, out of plane with the missile. 
The largest acceleration component comes from the elevator, 
perpendicular to the wings. By placing the wings parallel 
with the plane of the missile, the largest acceleration 
component is used to create the largest missile corrections, 
perpendicular to the plane of the missile. The plane of the 
missile is defined by three points: target position, 
missile position and the projected impact point. 

A graph of a target acceleration, while performing 
missile evasion, might look like Figure III-3. The pilot 
commands A(max sust) or A(max inst) for a short time, then 
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reduces the acceleration to zero to regain lost energy, 
before applying the acceleration again. This process may be 
repeated 2,3 or more times during the missile flight time. 
The graph attempts to incorporate transient response induced 
by system time delays, transient responses of the control 
surfaces and pilot tendencies for maneuvering and control of 
the target. 




Figure III-3 Target Acceleration 

The resultant average acceleration is non-zero, with a 
non-zero variance. The mean and variance of Figure III-3 
can be estimated by the parameters assumed in Figure III-2. 

Figure III-l and Figure III-2 show total aircraft 
acceleration. The parameters as seen by the missile, in 
antenna coordinates, will vary depending on the three 
dimensional maneuver employed. The missile tracking 

subsystem must be designed to handle the maximum 
acceleration possible for each of the orthogonal components 
of the reference frame. 

Proportional navigation missiles compensate for the 
applied target acceleration by decreasing the line of sight 
rate induced by the change of target velocities. With a 
good target model and detection of maneuvering effects, the 
missile guidance can predict target motion and position. 
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IV. 



SYSTEM MODEL 



From the section on the ideal missile, it was 
ascertained that target parameters are not constant if 
lateral acceleration is applied. This thesis will attempt 
to incorporate as much sensor information as is available in 
defining the system models for an optimum missile. Current 
radars allow the measurement of range, range rate and off 
boresight bearing error. Measurements taken by the radar 
are referenced to a radar axis system. In order for the 
missile to use the information supplied by the radar, a 
common reference frame must be established. 

A. COORDINATE SYSTEM 

Each entity in the missile-target intercept problem has 
its own coordinate system. The overall geometry as seen 
from an "eye in the sky" would view it in space coordinates. 
An observer on the ground would view it in earth 
coordinates. The launching aircraft, missile and target 
aircraft will view it in an individual coordinate system, 
referenced to that specific platform. Trying to equate each 
coordinate system is not an easy task but one which is done. 

By use of Euler angles any reference coordinate system can 
be related with Earth coordinates. By use of a directional 
cosine matrix transformation any reference coordinate system 
can be transformed to another reference coordinate system 
[Ref . 3] . 

The missile is concerned with flying to a point in space 
that will hopefully be occupied by the target at the 

completion of the intercept. The object is to guide the 

missile to the proper point where the target will be. The 
missile is concerned with its coordinate system and not 

that of the target. But on the missile itself there are 
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various coordinate system reference points. Each sensor has 
its own location on the missile and where it is mounted is 
its reference point. Any moving sensor, like the antenna, 
will have its special reference coordinate system. Missile 
parameters are normally referenced to the missile body frame 
of reference while target parameters are referenced to the 
antenna frame of reference. While very complicated, the 
frames of references can be transformed and equated . [Ref . 3] 
To simplify simulation and evaluation of desired 
parameters, an inertial frame of reference will be used 
which is centered at the radar antenna location. This 
simplification will aid in better evaluation of the effects 
of the target parameters and the missile guidance without 
encumbrance of transformation errors. Although the 

simplification assumes ideal missile parameters, time delays 
can be incorporated later to account for first order 
modeling of the missile. 



B. EQUATIONS OF MOTION 



In cartesian coordinates missile and target motion is 
described by the standard motion equation: 
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The equation is based on a fixed reference point. The 
orthogonal directions (Y and Z) will have the same equation. 

The antenna frame of reference uses polar coordinates 
which have the equations; 
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where r(t) = the radial component of motion 
e(t) = horizontal component of motion 
4>(t) = vertical component of motion 

When the reference point is not fixed, extra terms and 
cross coupling are introduced into equation dynamics. The 
Coriolis equation accounts for the moving reference point. 
Depicted in Figure IV-1 a change in the vector R is 
accounted for by both changes in the magnitude and the 
rotation effects by the moving reference point. Using the 
terms as defined in Figure IV-1 we can obtain the necessary 
equations to find r(t), 0{t) and <i>(t). For simplicity, only 
the derivation for r and © are shown with r and <i> 
relationships being a duality of derivation of r and ©. The 
simulation of Section VI will be two dimensional. 




The Coriolis equation to relate the time rate of change 
of the R vector to the rate of change in the r direction and 
the angular rate of motion is: 

R = r + wxR ( 11 ) 

where R = the directional vector 

r = the magnitude of the directional vector 
w = the angular rate of motion 
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The total change of the vector R is the sum of the change in 
the magnitude of R due to changes along the original vector 
R given by r and the angular rotation due to the moving 
coordinate frame, given by wxR. 

Utilizing the general rule for differentiating a vector^ 



an expression is obtained for the acceleration of the R 
vector . 

R = r + wxr + wxr + wxR (12) 
R = r + wxr + wxr + wxr + wxwxr (13) 
R = r + wxr + 2 (wxr) + wxwxr (14) 



where R = directional vector 

r = the acceleration of the magnitude of R 
wxr and wxr are perpendicular to the R vector 
wxwxr is centrifugal acceleration 



This equation gives the cross correlation of range and 
angle to implement in a second order model. Applying the 

rule of differentiating a vector again will yield the 
equations for a third order model. 

• • •• •• •• • 
d (R ) = r + w xr + wxr + 2 (wxr) + 2 (wxr ) + wxwxr 
^3 1 . ... 

+ wxwxr + wxwxr + wxR (15) 



d (R ) = r + w xr + wxr + 2 (wxr) 
dt . . ... 

+ wxwxr + wxwxr + wxr + wxwxr 
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+ 



2 (wxr ) + wxwxr 
wx2(wxr) + wxwxwxr 
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= r + 3 (wxr) + 3 (wxr ) + w xr + wxwxr 

+ 2 (wxwxr) + 3 (wxwxr) + wxwxwxr 



(17) 



where 



is the acceleration jerlc of vector R 



‘ rule for differentiating a vector 
+ wxA 

the total derivative is the sum of the time rate of change 
of the vector and the rotation of the vector. 
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r is the magnitude of R 

wxr is perpendicular to R 

wxwxr is in the negative direction of R 

wxwxwxr is perpendicular to R 

C. SECOND ORDER MODEL 



Beginning with equations 8-10, a second order state 
space system can be set up which would have the form: 
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(18) 

(19) 

( 20 ) 



The range portion of the second order system may be 
accomplished totally by the radar, since no other subsystems 
require the information. The radar receiver is designed to 
trac){ the target in the radial direction without the need 
for an additional filter. 

The more difficult state equations to implement in the 
missile are the angular directions. The second order, time 
invariant, constant velocity, zero acceleration, state 
feedbac)c model ma)ces the tracking much easier. The 
continuous system model can be given in a time derivative 
form as: 



X*+ ki-X + k2-X = 0 



(21) 



If the term given in the equation as X is actually the error 
of the angular position then 

X = (line of sight angle - antenna position angle) 
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where X can be directly measured by the antenna. The values 
of the k's in the time derivative equation depend on the 
designer and the response desired. For the simulations of 
the proportional navigation missiles of Section II and 
Section VI, ki = 20 and k 2 = 100 were used. These constants 
give a response time constant of 0.1 sec. 

The use of the coriolis equation to derive the second 
order model gives a time varying solution. Implementing 
time varying equations are difficult and often avoided by 
using a time invariant system and state feedback to cancel 
errors. The time variant space state model derived from 
equation 14 is shown below: 

R = r + wxr + 2(wxr) + wxwxr (22) 



Separating into orthogonal components of radial and 
transverse with scalar multiplication: 



aF = r + wwr (23) 

ar = wr + 2 (wr) (24 ) 

Rearranging into equations to implement into a system: 



r'= -02r + as (25) 

e'= -^r + ai (26) 

r r 

The state space model is difficult to represent unless 
divided into two channels with cross coupling given in the A 
matrix . 
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D. THIRD ORDER MODEL 

The time invariant model, derived from equations similar 
to those deriving the time invariant second order model 
(equations 18 and 19) , in state space form is given by: 
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Ranging may be accomplished by the radar receiver, as in 
the second order model, for similar reasons. Tracking 
angular positioning requires knowledge of the target angular 
acceleration as well as angular velocity. A filter is 
normally used to maintain a track of the target angular 
parameters. Common filters are a-(3 , Weiner, and Kalman. As 
compared in Reference 4, the Kalman is the best filter 
suited for air to air missiles, but also the most costly to 
implement. For the time invariant third order model a 
simple constant gain Kalman Filter can be used. The Kalman 
Filter will be discussed in the next section. 

The tine variant third order model is obtained from the 
second derivative of the coriolis equation derived in the 
previous section. 

R = r + 3{wxr) + 3{wxr) + w xr + wxwxr 

+ 2 (wxwxr) + 3 (wxwxr) + wxwxwxr (30) 

The R term is the change in the acceleration of the 
vector or a "jerk" term, a simple comprehension is the rate 
at which the pilot applies the commanded acceleration. 
Separating the equation into radial and tangential terms, 
the two orthogonal scalar equations are: 

3r = r - 3w^r - 3wwr (31) 
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3wr + 3wr + w r + r 



(32) 



at = 

Converting the equation into primary coordinate axis, 
the equations obtained are: 

*r‘= as + 302r + 3e(e‘)r (33) 

‘e‘= ai - 3(e')r - 38(r ’ ) - (34) 

r r 



The disassociated space state model looks like: 
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(35) 



(36) 



It is readily observed that a very high cross coupling 
of the radial and transverse components exists. The range, 
range rate and range acceleration are required to adequately 
compute the angular acceleration. The angular velocity and 
acceleration is required in computing the range 
acceleration. All of these quantities are time varying 
requiring a time varying filter to implement this model. 

The cubic term of angle rate in the angle channel is 
insignificant compared to the other terms and is neglected. 
The second order model uses the simplifying assumption of 
constant velocity and constant acceleration. For the full 
third order model, no simplifying assumptions will be made. 
This third order model should account for all of the cross 
coupling between the bearing and angle channels. 

A Kalman Filter can be employed to track the target in 
both range and bearing to implement the time variant third 
order model. The Kalman Filter will be discussed in the 
next section. 
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V. 



KALMAN FILTER 



Given a system model, where the plant can be modeled by 
a set of first order differential equations and the output 
can be measured, a set of state equations can be defined 
similar to; 

X=AX + BU + W (37) 

Y = C X + V (38) 

where X is the state vector 

Y is the system output vector 

A,B,C and D are matrices 

U is the system input 

W is plant disturbances 

V is measurement noise 

The system can be modeled in discrete time as: 

X{k+1) = ^ X(k) + r U(k) + W(k) (39) 

Y{k) = H X(k) + V{k) (40) 

A Kalman Filter is the best filter to track the output 
of a discrete system [Ref. 3]. The Kalman Filter equations 
are given as: 



X{“ |k ) 


= X{k |k- 


1 ) 


+ G {k ) • 1” 


^Y(k) - Y(“ |k- j ) 


ij 


X{X ^ » |k 


) = 4>-X{»‘ 1 


k ) 


+ r*u(k) 




(42) 


Y{»‘^ 1 |k 


) = H-X(“* 


‘ 1 


k ) 




(43) 


where 


X{k |k ) 


= 


the state 


estimate at time 


k given 






information 


through time k. 






X{‘'» 1 Ik ) 


= 


the state 


estimate at time 


k+1 given 



information through time k. 

Y{‘‘** |k) = the output estimate at time k+1 
given information through time k. 
G{k) = the filter gain at time k. 
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For linear, time invariant systems, the <{> and T matrices 

are easy to calculate and follow directly from the state 

At r At 

space model, where 4> = e and F = e dt . For non linear 
systems, an extended Kalman filter can be used. For the 
extended Kalman filter, the 4 > and F matrices are linearized 
about the projected operating point. One method of 

estimating the linearization is to take the partial 
derivatives of the non linear state space matrices; 



<I> 



F 



A 

X“ 

X = Xo 
U = Uo 
W = 0 



B 

U 



X = Xo 
U = Uo 
W = 0 



(44) 



(45) 



The gain matrix G(k) will vary with the parameters of 
the filter. G(k) is the weighting factor for the system 
error. The solution to the filter gain G(k) requires the 
solution of Riccati equations: 

G(k) = P (X |k-i ) HT 1^ H P (X |k- , ) HT + R(k) J ^ (46) 

P(X |k-i ) = 4> P(X |k-i ) 4>T + 6 Q 5T (47) 

P(X|k) =P(X|k-i) - G(k) H P(X|k-i) (48) 

where G(k) = Kalman Filter gain at time k 

P(X|k-i) = Covariance of predicted estimate 
R(k) = measurement covariance matrix, ECVV^} 

Q(k) = maneuver covariance matrix, ElUU^j 
P(x |k) = Covariance of filtered estimate 
5 = maneuvering weighting matrix 

A constant gain matrix can also be used in the Kalman 
filter. Instead solving the full Riccati equations for each 
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change of variables, a constant value is used throughout the 
problem. A constant gain matrix will simplify 
implementation of the Kalman filter. 

One implementation of the third order model, as 
discussed in the previous section, is to model the system as 
linear, time invariant, given by the space state model: 
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The Kalman Filter equations for this third order model 



are : 
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Using a Kalman Filter on this third order model is very 
simple and requires few on line calculations. The gain 
matrix can be considered either constant or time varying. 
If time varying gains are used, they can be computed off- 
line and stored in memory. The Filter then utilizes the 
precomputed gain schedule and can select a gain depending on 
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the accuracy of the filter at that time. If a maneuver 
causes the filter to loose accuracy then a higher gain term 
can be utilized. If constant gains are used then they must 
be high enough to compensate for any maneuver the target 
might make. A high gain matrix will make the missile more 
responsive to any unwanted noise terms in the system since 
the missile cannot distinguish a noise input from a target 
maneuver . 

As discussed previously, a Kalman Filter is not required 
for the range channel. The radar can measure range and 
range rate directly. Since the actual values of the range 
channel are not used by any other elements of the guidance 
subsystem, the radar is able to maintain its own tracking of 
the target in the center of the range gate, which has no 
consequences on the rest of the missile guidance. Some 
noise information can be gained when estimating the range 
channel with a Kalman Filter. A Kalman Filter is used for 
the range channel in the simulation of Section VI for 
completeness of simulation and practical experience. 

A second implementation of the third order model is by 
using the equations obtained through the coriolis equations. 
The disassociated state space model is given as: 
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(53) 



(54) 



Two possible ways to implement the Kalman filter are to 
create an extended Kalman filter by linearizing the $ 
matrix, reducing the time dependence and cross coupling of 
the range and bearing channel, or keep the cross coupling 
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components and have the Kalman Filter maintain the values of 
the time varying <t'. If is linearized, some of the cross 
coupling and time dependence lost by linearization can be 
compensated for by the maneuver covariance matrix Q. Given 
the target acceleration probability model. Figure III-l, the 
Q matrix can be calculated, as derived in Reference 3, as 
time varying and relates the cross coupling of the bearing 
and range channel as ; 
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where = acceleration variance 

As discussed in the reference the Q(3,3) element can be 
increased to make the missile gain matrix put more weight on 
any target acceleration elements. 

If the •t’ matrix is maintained as time varying and 
nonlinear then the Q matrix can be constant. The constant Q 
matrix can be calculated as: 

Q = K I (56) 

j 

where K = matrix gain 

I = identity matrix 

Since the reference deals extensively with the time 
invariant and the time varying Q matrix, this thesis will 
deal with the time varying and constant Q. 

If Figure III-4 is used to define the maneuver 
probability then a non-zero mean is established. The Q 
matrix maintains the same properties just discussed with a 
different calculation for o . The non-zero mean can be 
implemented by increasing Q(3,3), increasing the weighting 
matrix 5, or by not assuming the U(k) term is zero. 
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Kalman filter equations for this third order model are: 



xc* |k ) 
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X(“ |k-i ) 


+ G(k ) • 1^ Y(k) - Y(X |k_i ) j 


(57) 
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- G(k) -H-P(x |k_i ) 


(62) 
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VI . SIMULATION 



To aid in the efforts of simulation, the Dynamic 
Simulation Language (DSL) was used to integrate the 
equations of motion for the target and missiles, as well as 
the antenna angle channel. Two programs were written and 
are listed in Appendix B. The first program is the time 
varying third order model. The second program is the time 
invariant third order model and the second order model. 

The output of the simulation is a set of graphs to 
compare the three missile models and their effectiveness in 
tracking the target. 

The Kalman Filter is implemented in a Fortran subroutine 
at the end of the DSL main program. The basic filter 
equations used were described in the previous section. 

A. ASSUMPTIONS 

The following assumptions are made to simplify the 
implementation of the Kalman Filter and determine the 
effects of the time varying third order model. 

- "Ideal" missile autopilot. 

- Inertial cartesian reference frame for angle 
measurements . 

- Final portion of intercept only. 

- Cross coupled effects of missile motion on antenna 
stabilization system disregarded. 

- Missile initialized to collision course. 

- Missile constant speed of 1500 kts (Mach 1.5) or 
2500 ft/sec. 

- Target constant speed of 400 kts (Mach .75) or 667 
ft/sec . 

- Target lateral acceleration applied perpendicular 
to target velocity vector. 

- Angle of Attack not accounted for in velocity 
vector calculations. 

- Missile located at the center of the cartesian 
reference frame. 
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- Target located to the right of the origin of the 
reference frame. 

- No noise. 

B. INITIALIZATION 

The user is asked at the beginning of the program to 
establish the geometry by giving the target initial 
position, heading, speed and acceleration. The target 
heading is oriented relative to a vertical line, parallel to 
the y axis, defining the North or 000 heading. The initial 
missile headings are calculated for constant velocity, zero 
acceleration collision course with a time to go estimate of 
range/missile velocity. Antenna parameters are initialized 
to initial line of sight and zero angular rate. 

C. SECOND ORDER MODEL 

The second order model is implemented using a 
proportional navigation constant of four and two s-plane 
poles at s=-10. This gives the equation for angular 
acceleration as: 

p'= -20 B + 100 (LOS - (5) (63) 

where P = the antenna angular acceleration 
P = the antenna angle rate 
P = the antenna angle position 
LOS = the actual angle to the target 

D. THIRD ORDER MODEL 

The Kalman filter is used to implement the third order 
model. The DSL main program calls the Kalman Filter 
subroutine at a sampling time of h=0.01 seconds. The Kalman 
Filter is executed, then control is passed back to the DSL 
program. A proportional navigation constant of four is used 
as discussed in Section II. 
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1 . Time Invariant, Constant Phi Model 

The discrete time invariant third order model 
divided into two Kalman filters for range and bearing is: 



RNGC< ^ 1 |k ) = 


RPHI •RNG(“ |k ) 


(64) 


RNG (“ |k ) = RNG(“ 


|k - 1 ) + GR(k ) • [ DELR ] 


(65) 


S * 1 |k ) = SPHI • 


S (*' |k ) 


(66) 


S (“ |k ) = S (k jk- 1 


) + GS (k ) • SDEL ] 


(67) 



where 



RPHI = SPHI = 



1 .01 .0005 

0 1 .01 

0 0 1 



DELR = 



RM - RKPl 
RDM - RDKPl 



SDEL 



LOS - SKPl 
LOSD - SDKPl 



Using Matlab functions of Aker and Place, with 
eigenvalues of 0.5, 0.5 and 0.5 constant gain matrices were 
obtained for range and bearing, given by GR and GS , 
respectively . 



GR 



0.5 0.0125 

0.0025 1.0 
0.125 24.9 



GS = 



1 . 5 
12.5 
1250.0 



2 . Time Variant, Variable Phi Model 

The discrete, time variant third order model divided 
into two Kalman filters for range and bearing is: 

RNG(k |k ) = RNGC< |k - 1 ) + GR{k ) • [ DELR ] (68) 

RNG(k^i|k) = RPHI •RNG(*‘ |k ) 
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( 69 ) 



( 70 ) 



PR(“|k-i) = RPHI-PRC' |k- 1 ) -RPHIT + QR 

GR(k) = PR{“ |k-i ) -HRT r|HR-PR(x |k- 1 ) -HRT +RMCOV{kl 

PR(“|k) = PR{“|k-i) - GR(k) .HR-PR{“ |k- 1 ) 

S{X|k) = S(“|k-i) + GS(k)*[^ SDEL j 
S{“^ 1 |k ) = SPHI-S{“ |k ) 

PS{X|k-i) = SPHI-PS{“ |k-i ) -SPHIT + QS 



(73) 

(74) 

(75) 



GS(k) = PSC< |k-i ) *HST • I^HS'PSC^ |k-i ) *HST+RSCOV{k) j 

( 16 ) 

PS{k|k) = PSC^jk-i) - GS{k) •HS-PS{“ |k-i ) ( 11 ) 



The time variant model uses two different <t> 
matrices, one for range and the other for bearing. Two 

other matrices must be specified for each filter, the 
initial error covariance matrix and the target maneuvering 
covariance matrix. The two 4> matrices are: 



RPHI 



1 .01 .00005 

.00005*A 1+.00005*B .01 

.01*A .01*B 1+.00005*B 



(78) 



SPHI 



1 . 01 
0 1+.00005*C 

0 .01*C+. 00005*0 



.00005 

.01+. 00005*0 
1+. 01*0+. 00005* (C+0) 



(79) 



where 



A = -30 (e‘ ) 



B = -302 



C = -3r‘/r 
0 = -3r/ r 



The initial error covariance matrices are given as: 



PR(° |o ) 



500 0 0 

0 500 0 

0 0 500 
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PS (0 |o ) 



le4 

0 

0 



0 

le4 

0 



0 

0 

le4 



The maneuvering covariance matrix can either be time 
varying or constant as discussed previously. The 
maneuvering covariance matrix accounts for the capabilities 
of the target as discussed in Section IV. Reference 3 gives 
the derivation for the time varying matrix for a constant $ 
matrix. The <i> matrix also accounts for any time variance of 
the target model so the maneuvering covariance matrix can be 
constant. A constant matrix is assumed since <t> is time 
varying. The maneuvering covariance matrix (QR AND QS ) are 
given as: 



QR 



500 0 0 

0 500 0 

0 0 500 



QS 



.01 0 0 
0 .01 0 
0 0 .01 



The simulation calculates the range model then the 
bearing model. The results of each filter are used in the 
other filter to calculate the '!> values. 
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E. RESULTS 

Similar simulations of the ideal missile cases. Section 
II, were used to evaluate the missiles. The simulation 
consists of head-on, tail and beam aspects with OG, 3G and 
6G target acceleration. The results of the computer simu- 
lation is shown in graph form in Figure VI-1 through VI-31. 
The computer program listings are contained in Appendix B. 

1 . Gains 

Gain comparison plots are given in Figure VI-1 
through Figure VI-4. The resultant gains from the time 
variant, varying Phi third order missile are the same for 
each scenario. 

In predicting |k), the varying Phi model 
weights the range error by .5 and the range rate error by 0. 
The constant Phi model uses weights of .5 and .0125, 
approximately the same. Figure VI-1. 

Gains for predicting RD(‘‘*^ |u ) , for varying Phi, are 
.001 and .6 while those of the constant Phi are .0025 and 
1.0, Figure VI-2. Little emphasis is placed on the range 
error, because the radar is measuring range rate, with a 
higher weighting factor. 

In predicting RDD(*^*i |u ) , small gains are calculated 
by the varying Phi while the constant Phi model places a 
high emphasis on range rate error. The noise of the system 
will be noticed more in the prediction of RDD(‘‘** |k ) than 
the other parameters, because of the higher weighting 
factor . 

Bearing channel gains give unusual curves for the 
varying Phi model. Figure VI-4. There is little weight 
placed on non-observed parameters as in SG2 and SG3 , during 
most of the intercept, except in the initial stage and final 
stage. The gains are highest during the critical stages of 
the intercept. SGI is a constant 1.0 giving equal weighting 
to the current estimate and the error. The constant Phi 
model has higher gains giving more weight to any errors. 
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Figure Vl-1 Range Covariance Gains RGll and RG12 
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Figure VI-2 



Range Covariance Gains RG21 and RG22 
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Figure VI-3 Range Covariance Gains RG31 and RG32 
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Figure VI-4 
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Fioure VI-6 Line of Sicht Angle Plot for Head-on 
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Figure VI-7 Missile Commanded Acceleration Plot for Head-on 

Aspect No Target Turn 
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Fiaure VI-9 Line of Siant Angle Plot for Head-on 
Aspect 3G Target Turn 
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Figure VI-10 Missile Commanded Acceleration Plot for Head on 

Aspect 3G Target Turn 
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Fiaure VI-13 Missile Commanded Acceleration Plot for Head on 

Aspect 3G Target Turn 
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Figure VI-15 
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Figure VI-16 
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Figure VI-17 Position Plot for Tail 
Aspect 3G Target Turn 
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Figure VI-18 Line of Sight Angle Plot for Tail 
Aspect 3G Target Turn 
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Figure VI-IS Missile Commanded Acceleration Plot for Tail 

Aspect 3G Target Turn 
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VI -20 Position Plot for 
Aspect 6G Target Turn 
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REPRODUCED AT GOVERNMENT EXPENSE 
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Figure VI-21 



Line of Sight Anale Plot 
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Figure VI-22 Missile Corr.manded Acceleration Plot for Tail 

Aspect 6G Target Turn 
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REPRODUCED AT GOVERNMENT EXPENSE 
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Figure VI-23 Position Plot for Beain 
Aspect No Target Turn 
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Figure VI-24 Line of Sight Angle Plot for Beam 
Aspect IJo Target Turn 
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REPRODUCED AT GOVERNMENT EXPENSE 
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Figure VI-25 



Missile Cojrr.manded Acceleration Plot for Beam 
Aspect No Target Turn 
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Figure 



VI -26 Position Plot for 
Aspect 3G Target Turn 



Bean 
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Fiqure VI-27 Line of Siaht Anqle Plot for Beair. 
Aspect 3G Target Turn 
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Figure VI-28 Missile Commanded Acceleration Plot for Beam 

Aspect 3G Target Turn 
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REPRODUCED AT GOVERNMENT EXPENSE 












Fiaure VI-29 Position Plot for Beam 
Aspect dG Target Turn 
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Figure VI-30 Line of Siaht Angle Plot for Beam 
Aspect 6G Target Turn 
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REPRODUCED AT GOVERNMENT EXPENSE . 
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Figure VI-31 



Missile Commanded Acceleration Plot for Beam 
Aspect dG Target Turn 
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2 . 



Head-on Aspect 



Figures VI-5 through VI-16 are the results of 
missile guidance comparisons for head- on aspect initial 
condition with OG, 3G and 6G constant target acceleration. 
The missile begins at the origin of the graph. The target 
initial position is at x = 10000 ft, y = 1000 ft. Applied 
lateral target acceleration is away from the missile. 

a. OG Target Acceleration 

Flight paths for the three models are shown in 
Figure VI-5. All three paths appear to be the same, within 
the accuracy of the plotter. The line of sight angle. 
Figure VI-6 remains relatively constant for all three 
models, at the initial value, giving a constant bearing 
decreasing range, no manuever intercept. 

Commanded Missile Acceleration, Figure VI-7, 
shows that the second order model pulls more lateral 
acceleration than third order models. Higher gain terms 
make the constant Phi third order model erratic when 
compensating for initialization error. 

b. 3G Target Acceleration 

In the position plot. Figure VI-8, the second 
order model begins to lag the third order model . The 
lagging means the missile is slower to compensate for line 
of sight rates. This is further illustrated by Figure VI-9, 
the line of sight angle plot. The change in line of sight 
is greater for the second order model. 

Commanded acceleration. Figure VI-10 shows a 
larger increase for the second order model. The second 
order model requires approxima^'ely 7 . 5G to intercept a 3G 
target while the third order model require 4.66G. 
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c. 6G Target Acceleration 

Increased target acceleration increases the lag 
of the second order model. Figure VI-11. Magnitude of line 
of sight is larger for the second order model. Figure VI-12, 
but the missile is compensating for the errors. The two 
third order models are fairly close. The difference is only 
fractions of radians. 

The constant Phi third order model has less 
initial oscillations in commanded acceleration when the 
target acceleration increases. The second order model 
requires approximately IIG for a 6G target, 7 . 5G is required 
for the third order models. 

3 . Tail Aspect 

Figures VI-13 through Figure VI-21 are the results 
of missile guidance comiparisons for tail aspect initial 
condition with OG, 3G and 6G constant target acceleration. 
The missile begins at the origin of the graph. The target 
initial position is at x = 15000 ft, y = -500 ft. Applied 
lateral target acceleration is into the missile. 

a. OG Target Acceleration 

As shown in Figure VI-13, there is little 
difference in flight paths for the three missiles in the no 
target acceleration case. Line of sight angles remain 
constant. Figure VI-14, throughout the intercept. 

There is some slight difference in commanded 
accelerations. Figure VI-15, with the second order model 
being the smaller of the three models. 

b. 3G Target Acceleration 

As the target applies acceleration, the lag of 
the second order missile gives a shorter flight path, than 
the third order missiles. Figure VI-16. Line of sight 
angles are smaller for the second order model. Figure VI-17. 

Acceleration required for the second order model 
is 3 . 4G and 5.3G for the third order models. Figure VI-18. 
By lagging the other missiles, the second order model allows 
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the target to complete part of the intercept, which allows 
the missile to pull less acceleration. 

c. 6G Target Acceleration 

With higher target acceleration, the flight 
paths. Figure VI-19 have approximately the same differences 
as the 3G case. The second order model maintains a better 
flight path throughout the intercept. The slope of the line 
of sight curves. Figure VI-20, are higher with the second 
order model maintaining a smaller angle difference. 

Commanded acceleration is much smaller for the 
second order model than the third order models. Figure VI- 
21. To intercept the 6G target, the second order model 
requires 6 . 2G and the third order models require 8.1G. 

4 . Beam Aspect 

Figures VI-22 through VI-30 are the results of 
missile guidance comparisons for head- on aspect initial 
condition with OG, 3G and 6G constant target acceleration. 
The missile begins at the origin of the graph. The target 
initial position is at x = 15000 ft, y = 0 ft. Applied 
lateral target acceleration is away from the missile. 

a. OG Target Acceleration 

Figure VI-22 through VI-24, show the three 
missiles are practically the same for a no target 
acceleration, beam .aspect intercept. All three missiles 
maintain a constant bearing decreasing range, small 
acceleration intercept. 

b. 3G Target Acceleration 

Only slight differences are noticed when target 
acceleration is applied. Flight paths. Figure VI-25, shows 
very little deviation. The line of sight angle difference 
of .01 rad between the models is approximately . 575 . 
Acceleration is higher for the second order missile to allow 
it to fly the same path as the third order models. Figure 
VI-27. 
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c. 6G Target Acceleration 

Flight paths have a pronounced difference with a 
higher target acceleration, Figure VI-28. Line of sight 
angles, Figure VI-29, have larger magnitudes and higher 
slopes. Commanded acceleration increases dramatically, now 
25G is required of the second order model and 9.3G for the 
third order models . 



110 



VII , 



CONCLUSIONS 



In analyzing second and third order missile tracking and 

guidance subsystems, the following conclusions are made: 

Proportional navigation guidance is the optimum 
method for missiles, given current design 
tradeoffs . 

Target modeling is very difficult and requires 
the analysis or many factors. Acceleration 
probabilities make modeling easier, but the 
proper acceleration model must be chosen. 

Cross coupling between coordinate reference axis 
components does exist and gives errors if not 
accounted for in the system model. 

Kalman filters are the best predictors for 
airborne missiles, if one is required. 

Complete time varying third order models give 
better results than approximated linear, time 
invariant third order models. 

Only small differences are noticed in parameter 
values between second and third order models. 
Higher accelerations are required for the second 
order model. 

Second order missiles are better than third order 
missiles in tail aspect, constant acceleration 
intercepts . 

Implementation of a Kalman filter requires 
considerable amounts of computer resources, with 
limited time to complete the calculations. 

Some parameter terms are of the approximate order 
as system noise or non significant calculations. 

Some recommendations for future study and consideration 



are : 



A study of miss distance analysis for second and 
third order models. 

Analysis of the effects of the Q and P matrix 
initialization. 

- Analyze the target acceleration probability model 
to find optimum values to assign to the 
probability model. 

Determine missile cross couple effect of heading 
changes and autopilot torques on the sensor 
subsystem. 

Add noise to the system to determine the effects 
of the high gain terms on a noisy system. 



Ill 
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APPENDIX A. 



IDEAL MISSILE PROGRAM LISTING 



C MISSILE PROGRAM FOR FLIGHT PATH COMPARISON IN THE 

C THESIS, THIS PROGRAM IS DESIGNED TO FOLLOW THREE IDEAL 

C MISSILES USING DIFFERENT GUIDANCE TECHNIQUES TO TARGET 

C INTERCEPT. THE GUIDANCE TECHNIQUES ARE: PROPORTIONAL 

C NAVIGATION, PURE PURSUIT AND LEAD PURSUIT. PARAMETERS 

C ARE ALSO OBTAINED FOR A FOURTH MISSILE USING A DIRECT 

C FLIGHT PATH. 

INITIAL 

CONST 

G=32.2,D2R=.0175,K2F=1.66667,PITCH=2.7,PI=3.14159,RM1=20000 

K=0 

NN=0 

MISSX0=0.0 
MISSY0=0 . 0 
VM = 2500.0 

READ (2,5) VT, AT,THDG,TGTXO ,TGTY0 
5 FORMAT (F6.1,2X,F5.1,2X,F6.1,2(2X,F10.2) ) 

TGTV=VT*K2F 
TGTA=-AT*G 
THDG=THDG*D2R 
TGTVXO=TGTV*SIN (THDG) 

TGTVYO=TGTV*COS (THDG) 

C 

C INITIALIZATION OF DIRECT INTERCEPT MISSILE 

DXT=MISSXO 
DYT=MISSYO 
C 

C INITIALIZATION OF PROP NAV MISSILE CONSTANT VELOCITY, 

ZERO ACCEL 

R0=((TGTX0 - MISSX0)**2 + (TGTYO - MISSY0)**2 )**.5 
TTG0= RO/VM 

LOS = ATAN2 (TGTY0-MISSY0,TGTX0-MISSX0) 

PHDG=ATAN2 (TGTY0+TGTVY0*TTG0-MISSY0 , - 
TGTX0+TGTVX0*TTG0-MISSX0 ) 

BDO = TGTV*COS (THDG) / (COS (LOS) * RO ) 

BO = LOS 



C 

C INITIALIZATION OF LEAD PURSUIT MISSILE 

GO = BO 
GDO = 0 
LPHDGO = LOS 



Q -k -k 

METHOD RKSFX 

DERIVATIVE 

C 



c 

c 



c 

c 



c 

c 



TARGET POSITION UPDATING 
TGTHDG=INTGRL (THDG, (-1*AT) *PITCH*D2R) 
TGTAX=TGTA*COS (TGTHDG) 

TGTAY=-TGTA*SIN (TGTHDG) 

TVELX=INTGRL (TGTVXO ,TGTAX) 
TVELY=INTGRL (TGTVYO , TGTA^') 

XT=INTGRL (TGTXO , TVELX ) 

YT=INTGRL (TGTYO , TVELY) 

TGTHDG = AT AN2 (TVELX, TVELY) 

PROP NAV MISSILE POSITION UPDATING 
BDDOT = -20*BDOT + 100 * ( PNLOS-B ) 

BDOT = INTGRL (BDO , BDDOT) 

B = INTGRL (BO, BDOT) 

PNHDG = INTGRL (PHDG, 4*BDOT) 
PXM=INTGRL(MISSXO,VM*COS (PNHDG) ) 
PYM=INTGRL (MISSYO , VM* SIN (PNHDG ) ) 



PURE PURSUIT MISSILE 
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PPHDG = ATAN2 (YT-PPYM,XT-PPXM) 

PFXM = INTGRL (MISSXO,VM*COS (PPHDG) ) 

PPYM = INTGRL (MISSYO ,VM*SIN (PPHDG) ) 

C 

C LEAD PURSUIT MISSILE 

LPHDG 

ATAN2 (YT+TVELY*0. 5*TTG-LPYM,XT+TVELX*0. 5*TTG-LPXM) 

LPXM = INTGRL (MISSXO,VM*COS (LPHDG) ) 

LPYM = INTGRL(MISSY0,VM*SIN(LPHDG) ) 

C 

Q ******************* 

c 

DYNAMIC 

C 

C PROP NAV MISSILE GEOMETRY UPDATE 

PNAM = 4*BDOT*PNRD 

PNR= ( (XT-PXM) **2 + (YT-PYM) **2) ** . 5 
PNLOS = ATAN2 (YT-PYM, XT-PXM) 

PNRD=TVELX*COS (PNLOS) +VM* COS ( PNHDG-PNLOS ) 

PNLOSD = (-TVELX*SIN(PNLOS) -VM*SIN (PNHDG-PNLOS) ) /PNR 
C 

C PURE PURSUIT GEOMETRY UPDATE 

PPR = ( (XT-PPXM) * *2 + (YT-PPYM) **2) ** . 5 
PPLOS = ATAN2 (YT-PPYM, XT-PPXM) 

PPRD = TVELX*COS (PPLOS) - VM* COS ( PPHDG-PPLOS ) 

PPLOSD = -TVELX*SIN (PPLOS) /PPR 
C 

C LEAD PURSUIT GEOMETRY UPDATE 

LPR = ( (XT-LPXM) **2 + (YT- LPYM ) * * 2 ) * * . 5 
LPLOS = ATAN2 (YT-LPYM, XT-LPXM) 

LPRD = TVELX*COS (LPLOS) - VM* COS ( LPHDG-LPLOS ) 

LPLOSD = (-TVELX*SIN(LPLOS) 

VM*SIN(LPHDG-LPLOS) ) /LPR 
TTG =- LPR/LPRD 
C 

C DIRECT INTERCEPT MISSILE GEOMETRY UPDATE 

RD= ( (YT-MISSYO) **2 + (XT-MISSXO ) * * 2 ) * * . 5 
DXT , DYT , FON= CHECK ( RD , TIME , VM , XT , YT ) 

IF (PXM .GT. XT) CALL ENDRUN 
C 

IF (K .LE. 0.0) THEN 
NN=NN+1 

WRITE (31,50) XT, YT, PXM,PYM 
WRITE 32,50) PPXM , PPYM , LPXM , LPYM 
WRITE (33,51) TIME, PNLOS, PPLOS, LPLOS 
WRITE 34,52 TIME , PNLOSD , PPLOSD , LPLOSD 
WRITE (36,54) TIME , PKHDG , PPHDG , LPHDG 

50 FORMAT(4 (F10.2,2X) ) 

51 FORMAT (F5. 2, 3 (2X, F10.5) ) 

52 FORMAT ( F5 . 2 , 3 ( 2X , FIO . 5 ) ) 

54 FORMAT(F5.3,2X,3(F10.6,2X) ) 

K=10 

C 

C IF (PNR.LT.0.1*R0) THEN 

C K=1 

C ENDIF 

ENDIF 
K=K-1 
C 
C 

SAMPLE 

C SAVE (A) 0. 1, XT, YT, PXM, PYM, PPXM, PPYM, LPXM, LPYM 

C PRINT 1. 0, XT, YT, LPXM, LPYM, GDOT, LPLOS, LPHDG 
CONTROL FINTIM= 8.0,DELT=.01 
TERMINAL 

DHDG = ATAN2 (DYT-MISSY0,DXT-MISSX0) 

C GRAPH (A/A,DE=TEK618) XT 

(SC=1600,LO=0.00) ,YT(SC=750,LO=0.0,PO=16000) 

C GRAPH (A/A,OV) 

PXM ( SC=1600 , LO=0 . 0 , AX=OMIT) , PYM ( SC=750 , LO=0 ) 
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C GRAPH (A/A,OV) 

PPXM ( SC=1600 , LO=0 . 0 , AX=OMIT) , PPYM { SC=750 , LO=0 , AX=OMIT) 

C GRAPH (A/A,OV) 

LPXM ( SC=1600 , LO=0 . 0 , AX=OMIT) , LPYM ( SC=750 , LO=0 , AX=OMIT ) 

WRITE (2,15) DHDG,FON 

15 FORMAT (F10.7,2X,F5.2) 

WRITE (1,16) NN 

16 FORMAT (F4.0) 

END 

STOP 

FORTRAN 

SUBROUTINE CHECK (RD , TIME , VM , XT , YT , DXT , DYT , FON) 

REAL* 8 RD , TIME , VM , XT , YT , DXT , DYT , FON , DMISS 
DMISS = VM*TIME 
IF (DMISS .LT. RD) THEN 
DXT = XT 
DYT = YT 
FON = TIME 
END IF 
RETURN 
END 
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C THIS IS A DSL PROGRAM TO FIND THE LOS AND RANGE OF A 

C TARGET IN A CONSTANT G TURN FROM A MISSILE ON A 

C DIRECTPATH. 

C 

C 

INITIAL 

CONST 

G=32.2,D2R=.0175,PITCH=2.7 ,K2F=1 .66667 ,I=1.0,PI=3.14159,K=0 
READ (2,10) VT, AT,THDG,TGTXO,TGTYO 
READ (2,11) MISHDG,DONE 

10 FORMAT (F6.1,2X,F5.2,2X,F6.2,2(2X,F10.2) ) 

11 FORMAT (F10.7,2X,F5.2) 

C MISSILE PARAMETERS 

MISSXO = 0.0 
MISSYO =0.0 
VM = 2500.0 

C TARGET PARAMETERS 

TGTHDG=THDG*D2R 
TGTV=VT*K2F 
TGTG =-AT*G 

TGTTVELXO=TGTV* SIN (TGTHDG ) 

TGTTVELYO=TGTV*COS (TGTHDG) 

C EQUATIONS OF MOTION 

METHOD RKSFX 

DERIVATIVE 

TGTHDG=INTGRL (THDG, -AT*PITCH*D2R) 

TGTAX=TGTG*COS (TGTHDG) 

TGTAY=-TGTG*SIN (TGTHDG) 

TVELX=INTGRL (TGTTVELXO , TGTAX) 

TVELY=INTGRL (TGTTVELYO ,TGTAY) 

XT=INTGRL (TGTXO , TVELX) 

YT=INTGRL (TGTYO , TVELY) 

C MISSILE POSITION UPDATE 

XK = INTGRL (MISSXO ,VM*COS (MISHDG) ) 

YM = INTGRL (MISSYO, VM*SIN (MISHDG) ) 

DYI.'AMIC 

R= ( -XT-XM) **2 + (YT-YM) **2) ** . 5 
LOS = ATAN2 (YT-YM, XT-XM) 

RD = TVELX*COS (LOS) - VM*COS (MISHDG-LOS ) 

LOS: = ( TVELY*COS (LOS) - VM* SIN (MISHDG-LOS )) /R 
IF ^K.EQ.10) THEN 
MM = MM+1 

WRITE (37,15) XM,YM 

WRITE (38,16) TIME, LOS, LOSD, MISHDG 

15 FORi.AT (2(F9.2,3X)) 

16 FORKAT (F5.2,2X,2(F10.5,2X) ,F7.4) 

K=C 

END IF 
K=p;-i 

IF (TIME .GT. DONE) CALL ENDRUN 

SAMPLE 



CON’TROL FI2;TIM=10.0,DELT=.01 
C PRINT 1.0,XT,YT,XM,YM,LOS,R 

C SAVE (D) 0.1,XT,YT,XM,YM,R,LOS 

TERMINAL 

READ (1,26) NN 
WRITE (1,26) MM 
26 FORMAT (F5.0) 

C GRAPH (A/D,DE=TEK618) XT ( SC=1 500 , LO=0 ) , YT ( PO=l 5000 ) 

C GRAPH A/D,DE=TEK618,OV) XM ( SC=1500 , LO=0 , AX=OMIT) , YM 

C GRAPH (A/D,DE=TEK618) TIME, LOS 

C GRAPH (A/D,DE=TEK618) TIME , R 

END 

STOP 
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APPENDIX B. THIRD ORDER SIMULATION PROGRAM LISTING 



C PROP NAV MISSILE PROGRAM FOR THESIS 

INITIAL 

D DIMENSION 

PS (3, 3) , PR (3, 3) ,RMCOV(2,2) , RNG ( 3 ) ,S(3) ,DELR(3) ,SG(3) 

D DIMENSION RG{3,3) 

K=0 

MM=0 

HH=0 



C 

METHOD RKSFX 

CONST G=32 . 2, D2R=. 0175, K2F=1 .66667 
TK = .01 
MISSX0=0.0 
MISSY0=0 . 0 
VM = 2500.0 
AMO = 0.0 

READ (2,10) VT,AT,THDG,TGTXO,TGTYO 
10 FORMAT (F6.1,2X,F5.1,2X,F6.1,2X,F10.2,2X,F10.2) 

TGTV=VT*K2F 
TGTA= (-1*AT) *G 
THDG=THDG*D2R 
TGTVXO=TGTV*SIN(THDG) 

TGTVYO=TGTV*COS (THDG) 



C 

C 



C 

c 



c 

c 



c 

c 

MATRIX 



INITIAL PS (0/0) MATRIX 
PS(l,l)=1.0E+4 
PS 1,2 =0.0 
PS(1,3)=0.0 
PS (2 , 1) =PS (1 , 2) 
PS(2,2)=1.0E+4 
PS (2 , 3) =0. ■ 

PS 3, 1 =PS 
PS 3 , 2 =PS 




INITIAL PR(0/0) MATRIX 



i; 

2; 



PR (1 
PR(1 
PR (1, 3) 
PR(2, 1) 
PR (2, 2) 
PR (2, 3) 
PR (3,1 
PR(3,2) 
PR (3, 3) 



500 

0 

0 

PR(1 
500 
0 . 0 
PR(1, 
PR(2, 
500 



INITILIZE THE RANGE MEASUREMENT COVARIANCE MATRIX 
RMCOV(l,l) =0.0 
RMCOV 1,2 =0.0 

RMCOV(2,l) = RMCOV (1,2) 

RMCOV (2, 2) =0.0 

INITIALIZE THE BEARING MEASUREMENT NOISE COVARIANCE 



SMCOV = 0.0 
C 

C INITIALIZATION OF PROP NAV MISSILE CONSTANT VELOCITY, 

ZERO ACCEL 

LOS = ATAN2(TGTY0 - MISSYO , TGTXO - MISSXO) 
R=((TGTX0 - MISSX0)**2 + (TGTYO - MISSY0)**2 )**.5 
TTG0= R/VM 

PHDG=ATAN2 (TGTY0+TGTVY0*TTG0-MISSY0,TGTX0+- 
TGTVX0*TTG0-MISSX0) 

VMXO = VM*COS(PHDG) 
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VMYO = VM*SIN(PHDG) 

RKFi = R 

RDKFl = -VM*COS (PHDG-LOS) + TGTV* SIN (THDG ) /COS ( LOS ) 

RDDKPl = 0 
SKPl = LOS 

SDKPl = (TGTVYO/COS (LOS) - VM* SIN ( PHDG-LOS )) /R 
SDDKPl = 0 
BO = LOS 
EDO = 0 
C 

RNG(l) = RKPl 
RNG^) = RDKPl 
RNG(3) = RDDKPl 
S(l) = SKPl 
S 2) = SDKPl 

S(3) = SDDKPl 

C 

DERIVATIVE 

C 

C TARGET POSITION UPDATING 

TGTHDG = ATAN2 (TVELX,TVELY) 

TGTAX = TGTA*COS (TGTHDG) 

TGTAY= (-l‘TGTA) *SIN (TGTHDG) 

TVELX = INTGRL(TGTVXO, TGTAX) 

TVELY=INTGRL (TGTVYO ,TGTAY) 

XT = INTGRL(TGTXO, TVELX) 

YT=INTGRL (TGTYO , TVELY) 

C 

C 

C THIRD ORDER PROP NAV MISSILE POSITION UPDATING 

BDOT = INTGRL(BDO,BDDOT) 

B = INTGRL (BO , BDOT) 

AM = U 

MVELX = INTGRL (VMXO, -AM*SIN (PNHDG) ) 

MVELY = INTGRL (VMYO, AM* COS ( PNHDG ) ) 

PKHDG = AT AN2 (MVELY, MVELX) 

PXM= INTGRL (MI SSXO , MVELX) 

FYM= INTGRL (MISSYO , MVELY) 

C 

DYNAMIC 

C 

C THIRD ORDER PROP NAV MISSILE GEOMETRY UPDATE 

RM = ((XT-PXM)**2 + (YT-PYM) **2) »* . 5 

LOS = ATAN2 (YT-PYM, XT-PXM) 

RDOTM = TGTV*SIN (TGTHDG) /COS (LOS) 

VM*COS (PNHDG-LOS) 

RDDOTM = TGTA*COS (TGTHDG-LOS) 

LOSD = ( TVELX*SIN (LOS) + VM* SIN ( PNHDG-LOS )) /RM 

LOSDB = TGTAY*COS (LOS) /RM 
C COMPUTE THE ERROR TERMS 

DELR(l) = RM - RKPl 
DELR(2) = RDOTM - RDKPl 
SDEL = LOS - SKPl 
C 
C 

CALL 

KALMAN (RNG,RM, RDOTM, RDDOTM, RK,RDK,RDDK, RKPl, RDKPl, RDDKPl, . . . 



K, DELR, PR, RMCOV,S, LOS, LOSD, LOSDD,SK, SDK, SDDK, TIME, SKPl, . . . 
SDKPl , SDDKPl , SDEL, SMCOV, PS ,TK,RG, SG,HH) 
BDDOT=SDDKP1+10* ( SDK-BDOT) +33 . 33333* (LOS-B) 

U = -4*BDOT*RDOTM 
C 

C SAVE VALUES OF THE GAIN MATRICES 

GR11=RG (1 , 1) 

GR12=RG 1 , 2 
GR21=RG (2 , 1) 

GR22=RG (2, 2) 
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GR31=RG(3,1) 

GR32=RG(3,2) 

GS1=SG (1) 

GS2=SG (2) 

GS3 = SG (3) 

IF { PXM .GT. XT ) THEN 
CALL ENDRUN 
ENDIF 

IF (K .LE. 0) THEN 

WRITE (40,20) XT , YT , PXM , PYM 

WRITE 42,21 TIME,GR11,GR12,GR21,GR22,GR31,GR32 

WRITE (44,22) TIME , GSl , GS2 , GS3 
WRITE (46,23) TIME , LOS , LOSD , LOSDD , U 
K=10 



20 
21 
22 
23 
SAMPLE 
C 



IF (RM .LT. .1*R) THEN 
K=3 
ENDIF 
ENDIF 
K= K“ 1 

FORMAT(4 (2X,F10.2) ) 
FORMAT (F5.2,6(1X,F10. 
FORMAT (F5.2,3(2X,E12. 
FORMAT (F5.2,3 (2X,E11. 



4) 

5) ) 

4) ,2X,E14.6) 



C 

C 

c 

c 

c 

c 

c 

c 

c 



STATEMENTS TO SAVE DATA FOR USE WITH GRAFAEL 
SAVE (A) 0.1,XT,YT,PXM,PYM 

(B) 0.1, LOS,B,SK 

(C) 0.1, LOSD, BDOT, SDK 
p) 0.1,SDEL 

(E) 0.1,RM,RK,RDDOTM,RDDK 
0.1,U,BDDOT 
0.1,GS1,GS2,GS3 
0.1,GR11.GR12,GR21,GR22 



SAVE 

SAVE 

SAVE 

SAVE 

SAVE 

SAVE 

SAVE 

PRINT 



if 



G ) 
(H) 



0 . 1 ' RM , RKP 1 , RK , RDOTM RDKPl , RDK , RDDOTM , RDDKP 1 , RDDK 



PRINT 0.1, PNHDG , LOS , B , SK , LOSD , BDOT , SDK , LOSDD , BDDOT , SDDK 



CONTROL FINTIM=10.0,DELT=.01 
TERMINAL 

WRITE (1,30) MM 
30 FORMAT (F6.1) 

C STATEMENTS FOR PLOTTING WITH GRAFAEL 

C GRAPH (A/A,DE=TEK618) 

(SC=1600,LO=0.0) ,YT(SC=500,PO=16000) 



XT 



C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

END 
STOP 
FORTRAN 



GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 



(A/A,OV) PXM (SC=1600 , AX=OMIT) . PYM(SC=500) 
(B/B,DE=TEK618) TIME , LOS ( SC= . 025 , LO=- . 1 ) 

B/B,OV) TIME(AX=OMIT) , B ( PO=7 . 5 , SC= . 025 , LO=- . 1 ) 
(B/B,OV) TIME(AX=OMIT) , SK ( AX=OMIT , SC= . 0^5 , LO=- . 1 ) 
(C/C,DE=TEK618) TIME, LOSD 
(C/C,OV) TIME(AX=OMIT) , BDOT ( PO=7 . 5 ) 

C/C,OV) TIME(AX=OMIT) , SDK ( AX=OMIT) 
(D/D,DE=TEK618) TIME,SDEL 

(E/E,DE=TEK618) TIME,RM(SC=2000 .0,LO=0 .0) 

(E/E,OV) TIME (AX=OMIT) ,RK(SC=2006.0,LO=0.0) 

F/E , DE=TEK61 8 ) TIME , RDDOTM ( SC=500 . 0 ) 

(F/E,OV) TIME(AX=OMIT) , RDDK ( AX=OMIT) 
G/F,DE=TEK618) TIME,U 
(H/G,DE=TEK618) TIME, GSl 
(I/G,DE=TEK618 TIME,GS2 
(J/G,DE=TEK618) TIME,GS3 
K/H,DE=TEK618 TIME,GR11 
(L/H,DE=TEK618) TIME,GR12 
M/H,DE=TEK618 TIME,GR21 
(N/H,DE=TEK618) TIME,GR22 
(0/F,DE=TEK618) TIME, BDDOT 
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SUBROUTINE KALMAN ( RNG , RM , RDOTM , RDDOTM , RK , RDK , RDDK , 

RKPl , RDKPl , RDDKPl , K , DELR, PR , RMCOV , S , LOS , LOSD , LOSDD , 

SK , SDK , SDDK , TIME , SKPl , SDKPl , SDDKPl , SDEL , SMCOV , 

PS ,TK,RG, SG,HH) 

C SUBROUTINE TO ITERATE A KALMAN FILTER FOR RANGE 

C VARIABLES 

C GIVEN THE COVARIANCE MATRIX AND OBSERVATIONS 

REAL* 8 RNG (3) , RM , RDOTM , RDDOTM , RK , RDK , RDDK , RKPl , 

* RDKPl , RDDKPl , DELR ( 3 ) , PR ( 3 , 3 ) , RMCOV ( 3 , 3 ) , S ( 3 ) , 

* LOS , LOSD , LOSDD , SK . SDK , SDDK , SKPl , SDKPl , SDDKPl , SDEL , 

* SMCOV, TIME, PS(3,3) ,TK,TKSQ,A,B,C,D, 

* TEMPI (3), 
TEMP2 (2,2) , TEMP3 (3,3), RPHI (3,3), COVR (3,3) , 

* DET,SCOV,SPHI(3,3) , RG ( 3 , 3 ) , SG ( 3 ) , QR ( 3 , 3 ) , QS ( 3 , 3 ) 

* TKSQ = TK*TK 



C 

C 

OF 

C 



FIND THE NEW VALUES OF RPHI FROM THE PREVIOUS VALUES 

SIGMA 

MATRIX 

A = -3*SDKP1*SDDKP1 
B = -3*SDKP1*SDKP1 



RPHI 
RPHI 
RPHI (1 , 3 
RPHI (2,1 
RPHI (2,2 
RPHI (2 
RPKI (3 
RPHI (3 
RPHI ( 3 



3) 

1 

2 

3 



C 

C FIND 

RPHI*PR (K-l/K-1) *RPHI 

CALL ZERO(TEMP3, 3) 



1 

TK 

. 5*TKSQ 
. 5*A*TKSQ 
1 + .5*B*TKSQ 
TK 

A*TK 

B*TK + .5*A*TKSQ 
1 + .5*B*TKSQ 

THE PROJECTED COVARIANCE PR(K/K-1) 



C 



102 

103 

104 
C 

c 



105 

106 
107 
C 

C 



C 



110 

111 



DO 104 L=l,3 
DO 103 M=l,3 
DO 102 N=l,3 

TEMP3(L,M)=TEMP3(L,M) + RPHI ( L , N) * PR (N , M ) 
CONTINUE 
CONTINUE 
CONTINUE 

CLEAR THE OLD COVARIANCE MATRIX 
CALL ZERO (PR, 3) 

MULTIPLY BY RPHI TRANSPOSE 
DO 107 L=l,3 
DO 106 M=l,3 
DO 105 N=l, 3 

PR(L,M)=PR (L,M)+TEMP3 (L,N) *RPHI (M,N) 
CONTINUE 
CONTINUE 
CONTINUE 

DEFINE THE Q MATRIX OF MANEUVER COVARIANCE 



QR(1,1) = 500 
‘ " ■ 0.0 
0.0 
0.0 
500 
0.0 

QR ( 1 , 3 ) 

, OR ( 2 , 3 ) 

3) = 500 

NOW ADD TO THE COVARIANCE MATRIX 
DO 111 L=l,3 
DO 110 M=l,3 

PR(L,M) = PR(L,M) + QR(L,M) 
CONTINUE 
CONTINUE 



QR 

QR 

QR 

QR 

QR 

QR 

QR 

QR 



§! 
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C NOW PR IS THE COVARIANCE OF THE PREDICTED ESTIMATE ' 

PR(K/K-1) 

C FIND THE ESTIMATE OF RANGE MATRIX AT STEP K 

C 

C ZERO A TEMPORARY MATRIX 

CALL ZERO (TEMP2,2) 

DO 121 L=l,2 
DO 120 M=l,2 

TEMP2(L,M) = PR(L,M) + RMCOV(L,M) 

120 CONTINUE 

121 CONTINUE 

DET = TEMP2(1,1) *TEMP2(2,2) - TEMP2 { 1 , 2) »TEMP2 { 2 , 1 ) 

COVRd.l) = TEMP2(2,2)/DET 
C0VR(1-2 = (-lj*TEMP2(l»2)/DET 

COVR(2,l) = (-1) *TEMP2(2,1)/DET 
COVR(2,2) = TEMP2(1,1)/DET 
C HERE COVR = (HPH + R) INVERSE 

C ZERO A TEMPORARY MATRIX 

CALL ZERO (RG,3) 

DO 132 L=l,3 
DO 131 M=l,2 
DO 130 N=l,2 

RG (L,M) =RG (L,M) + PR(L,N) * COVR(N,M) 

130 CONTINUE 

131 CONTINUE 

RG(L,3) = 0 

132 CONTINUE 

C RG = PH (HPH + R) (INVERSED) H 

C ZERO A TEMP MATRIX FOR THE RNG MATRIX 





DO 140 L=1.3 
TEMPI (L) =6 


140 


CONTINUE 



C 



DO 142 L=l,3 
DO 141 M=l,2 





TEMPKL) = TEMPKL) 


+ RG(L,M) 


* DELR(M) 


141 


CONTINUE 




142 


CONTINUE 
DO 143 N=l,3 








RNG(N) = TEMPKN) + 


RNG (N) 




143 


CONTINUE 






C 


SAVE THE VALUES OF RANGE 


MATRIX AT 


STEP K 



RK = RNG(l) 
RDK = RNG (2) 
RDDK = RNG (3) 



C 

C ZERO THE OLD RANGE TEMPORARY MATRIX 

DO 150 N = 1.3 
TEMPI (N) = 6 

150 CONTINUE 

C FIND THE ESTIMATE OF THE STEP K+1 FOR THE RANGE MATRIX 

DO 152 L = 1, 3 
DO 151 M = 1,3 

TEMPKL) = TEMPKL) + RPHKL.M) * RNG(M) 

151 CONTINUE 

152 CONTINUE 
C 

C SAVE THE VALUES OF RNG(K+1/K) 

DO 153 N=l,3 

RNG(N) = TEMPI (N) 

153 CONTINUE 
C 

RKPl = RNG(l) 

RDKPl = RNG (2) 

RDDKPl = RNG (3) 

C FIND THE COVARIANCE OF FILTERED ESTIMATE 

C RG = PH (HPH+R) INVERSED H 

C THEREFORE P(K/K) = P(K/K-1) + RG*P(K/K-1) 
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161 

162 

163 

C 



164 

165 
C 

C 

c 

c 

VAR 

C 

c 

c 



CALL ZER0(TEMP3, 3) 

DO 163 L=l,3 
DO 162 M=l,3 
DO 161 N=l,3 

TEMP3(L,M) = TEMP3(L,M) + RG ( L , N) * PR (N , M) 

CONTINUE 
CONTINUE 
CONTINUE 

DO 165 L=l,3 
DO 164 M=l,3 

PR(L,M) = PR(L,M) - TEMP3(L,M) 

CONTINUE 
CONTINUE 

NOW PR IS THE COVARIANCE OF FILTERED ESTIMATE P(K/K) 

SUBROUTINE TO ITERATE A KALMAN FILTER FOR SIGMA 
TABLES 

GIVEN THE COVARIANCE MATRIX, OBSERVATION. 

WHERE H= (1 0 0) 

CALCULATE THE NEW SPHI MATRIX 
C = -3*RDDK/RK 
D = -3*RDK/RK 



C 

C 

c 

c 

c 

c 



2C . 
2C ^ 

2c : 

c 



203 

204 

205 
C 



SPHI (1 
SPHI {1 
SPHI (1 
SPHI (2 
SPHI (2 
SPHI 2 
SPHI (3 
SPHI 3 
SPHI 3 



1 ) 

2 

3 

1 

2 

3 

1 

2 

3 



1 
TK 

. 5*TKSQ 
0 

1 + .5*C*TKSQ 
TK + .5*D*TKSQ 
0 

C*TK + .5*C*D*TKSQ 
1 + D*TK + . 5* (C+D) *TKSQ 



CALCULATE THE NEXT PROJECTED P MATRIX FOR NEXT STEP 

T 

PS (K^ K-1 ) =SPHI*PS*SPHI 

CLEAh A TEMPORARY MATRIX 
CALL ZERO(TEMP3 , 3) 

EC 202 L=l,3 
DO 201 M=l,3 
DC 200 N=l,3 

TEMP3 (L,M) =TEMP3 (L,M) + SPHI ( L , N ) * PS (N , M) 
CONTINUE 
CONTINUE 

c:ntinue 

ZERC THE OLD PS MATRIX 
Cr-.LL ZERO (PS, 3) 

DO 205 L=l,3 
DO 204 M=l,3 
DO 203 N=l,3 

PS (L,M)=PS (L,M)+TEMP3 (L,N) *SPHI (M,N) 

CONTINUE 

CONTINUE 

CONTINUE 

DEFINE THE Q MATRIX OF MANEUVER COVARINCE 
QS(1,1) = 0.01 

g£(l,2) = ‘ ‘ 



S(l,3 
_S (2,1 
OS (2,2) 
OS (2,3) 
QS 3,1 
OS (3,2) 



0.0 
0.0 
= OS (1,2) 
= 0.01 
= 0.0 
= QS(1 
= OS (2 
0 . 01 
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S(3,3) 

O 211 L=l,3 
DO 210 M=l,3 

PS(L,M) = PS(L,M) + QS(L,M) 
CONTINUE 
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211 

C 

C 

C 

C 

C 



C 



C 

C 



c 

c 

c 

c 

c 

c 



220 



221 

222 

C 

C 



223 

C 

C 



C 

C 

C 



450 

451 
C 



CONTINUE 

T 

CALCULATE (HPH + R) INVERSE 



SCOV = PS (1,1) + SMCOV 

SG(1) = PS(l,l)/SCOV 

SG(2) = PS 2,1 /SCOV 
SG(3) = PS{3,l)/SCOV 

NOW FIND THE CURRENT VALUES OF THE SIGMA MATRIX 
S(l) = S(l) + SG(1)*SDEL 
S (2) = S (2 + SG (2 *SDEL 

S(3) = S(3) + SGU)*SDEL 

STORE THE SIGMA MATRIX FOR USE IN THE PROGRAM 
SK = S(l) 

SDK = S(2) 

SDDK = S(3) 

FIND THE NEXT VALUES OF THE SIGMA MATRIX 

S(K+1) = SPHI * S(K) 

ZERO A TEMPORARY MATRIX 

DO 220 L=l,3 

TEMPI (L) = 0 
CONTINUE 
DO 222 L=l,3 

DO 221 M=l,3 

TEMPI (L) =TEMP1 (L) 

CONTINUE 
CONTINUE 



+ SPHI(L,M) * S(M) 



INPUT BACK INTO SIGMA MATRIX 
DO 223 N = 1,3 
S (N) = TEMPI (N) 

CONTINUE 

STORE THE VALUE OF THE S MATRIX 
SKPl = S(l) 

SDKPl = S(2) 

SDDKPl = S(3) 



NOW FIND THE P MATRIX AT STEP K 



IF 



PS (3, 3 
PS (3,2 
PS 3,1) 
PS(2,3)^ 
PS (2, 2) 
PS(2,1) 
PS 1,3 
PS(1,2)= 
PS(1,1)= 
(HH .LE, 
HH=50 




PS (1, 
0 . 0 ) 



*SG(3) 
*SG (3) 
*SG (3) 
*SG(2) 
*SG(2) 
*SG (2) 
»SG 1) 
*SG (1) 
*SG (1) 



WRITE (9,*) 'TIME IS ’,TIME 
WRITE 9,* ‘A,B,C,D,’ 

WRITE (9,*) A,B 
WRITE 9,* C,D 

WRITE (9,*) 'RANGE PHI MATRIX' 

DO 450 L=l,3 

WRITE (9,*) (RPHI(L,M), M=l,3) 
CONTINUE 

WRITE (9,*) 'SIGMA PHI MATRIX' 

DO 451 L=l,3 

WRITE (9,*) (SPHI(L,M), M=l,3) 
CONTINUE 

OUTPUT THE COVARIANCE MATRICES AT STEP K 
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no 



401 

402 

300 

301 



WRITE (9,*) 'THE RANGE COVARIANCE MATRIX IS:' 
DO 401 M=l,3 

WRITE (9,*) (PR(M,N) , N=l,3) 

CONTINUE 

WRITE (9,*) 'THE BEARING COVARIANCE MATRIX IS 
DO 402 M=l,3 

WRITE (9,*) (PS{M,N) , N=l,3) 

CONTINUE 
END IF 

HH=HH-1 

RETURN 

END 

SUBROUTINE ZERO (A, N) 

CLEAR A TEMPORARY MATRIX 

REAL* 8 A (3, 3) 

DO 301 L=1,N 

DO 300 M=1,N 
A (L,M) =0 
CONTINUE 
CONTINUE 
RETURN 
END 



124 



C THIRD ORDER MISSILE SIMULATION USING CONSTANT GAINS 

FOR 

C THE KALMAN FILTER. A SECOND ORDER PROP NAV REFERENCE 

C MODEL IS SIMULATED FOR PLOTTING AND PARAMETER 

COMPARISONS. 

INITIAL 

D DIMENSION 

RNG { 3 ) , S ( 3 ) , DELR ( 2 ) , GS ( 3 ) , GR ( 3 , 3 ) , RPHI (3,3), SPHI (3,3) 

K=0 

NN=0 



C 

C 

C 

METHOD RKSFX 

CONST G=32.2,D2R=.0175,K2F=1.66667 
TK = 0.01 
MISSX0=0.0 
MISSY0=0.0 
VM = 2500.0 
AMO = 0.0 

READ (2,10) VT, AT,THDG,TGTX0 ,TGTY0 
10 FORMAT (F6.1,2X,F5.1,2X,F6.1,2(2X,F10.2) ) 

TGTV=VT*K2F 
TGTA=-AT*G 
THDG=THDG*D2R 
TGTVX0=TGTV*SIN(THDG) 

TGTVY0=TGTV*COS (THDG) 



C 

c 



c 

c 



c 



c 

c 



INITIALIZE THE 
RPHI (1 , 1) 
RPHI (1,2) 
RPHI (1 , 3) 
RPHI 2,1 
RPHI (2, 2) 
RPHI (2, 3 
RPHI (3, 1) 
RPHI (3, 2 
RPHI (3, 3) 



RANGE PHI MATRIX 
= 1 
= TK 

= .5*TK*TK 
= 0 
= 1 
= TK 
= 0 
= 0 
= 1 



INITIALIZE THE 
SPHI (1,1) = 
SPHI(1,2) = 
SPHI 1,3) = 
SPHI(2,1) = 
SPHI(2,2) = 
SPHI(2,3) = 
SPHI 3,1) = 
SPHI(3,2) = 
SPHI(3,3) = 
CONSTANT STEADY 



li:i! : 



GR 
GR 
GR ( 1 , 3 ) 
GR(2,1) 
GR 2,2) 
GR ( 2 , 3 ) 
GR 3,1 
GR(3, 2) 
GR ( 3 , 3 ) 



BEARING PHI MATRIX 
1 

TK 

.5*TK*TK 

0 

1 

TK 

0 

0 

1 

STATE GAIN VALUES RANGE 



5 

0.0125 
0 
0 
1 
0 

0 . 

24 
0 



0025 

0 

1250 

.9991 



CONSTANT STEADY STATE GAIN VALUES BEARING 
GS(1) = 1.5 
GS 2 = 12.5 

GS(3) = 1250.0 



C 

C INITIALIZATION OF PROP NAV MISSILE CONSTANT VELOCITY, 

Z£RO ACCEL 

LOS = ATAN2(TGTY0 - MISSYO , TGTXO - MISSXO) 
R=((TGTX0 - MISSXO) **2 + (TGTYO - MISSYO) **2 )**.5 
TTG0= R/VM 



125 



PHDG=ATAK2 (TGTYO+TGTVYO *TTGO-MISSYO , TGTXO+TGTVXO ‘TTGO-MISS- 
XO) 

VMXO = VM*COS(PHDG) 

VMYO = VM*SIN(PHDG) 

RKPl = R 

RDKPl = -VM*COS (PHDG-LOS) + TGTV*SIN (THDG) /COS (LOS) 

RDDKPl = 0 
SKPl = LOS 

SDKPl = (TGTVYO/COS (LOS) - VM* SIN ( PHDG-LOS ) ) /R 
SDDKPl = 0 
BO = LOS 
BDO = 0 
C 

RNG(l) = RKPl 
RNG{2) = RDKPl 
RNG(3) = RDDKPl 
S(l) = SKPl 
S(2) = SDKPl 

S(3) = SDDKPl 

C 

C INITIALIZATION OF SECOND ORDER PROP NAV REFERENCE 

MODEL 

GO=LOS 
GDO = 0 
C 

DERIVATIVE 

C 

C TARGET POSITION UPDATING 

TGTHDG = ATAN2 (TVELX,TVELY) 

TGTAX=TGTA*COS (TGTHDG) 

TGTAY= (-1*TGTA) ‘SIN (TGTHDG) 

TVELX=INTGRL (TGTVXO , TGTAX) 

TVELY=INTGRL ( TGTVYO , TGTAY) 

XT=INTGRL ( TGTXO , TVELX ) 

YT=INTGRL (TGTYO , TVELY) 

C 

C 

C THIRD ORDER PROP NAV MISSILE POSITION UPDATING 

BDOT = INTGRL (BDO ,BDDOT) 

B = INTGRL (BO, BDOT) 

AM = U 

MVELX = INTGRL (VMXO , -XM‘SIN (PNHDG) ) 

MVELY = INTGRL ( VMYO, AM»COS (PNHDG) ) 

PNHDG = AT AN2 (MVELY , MVELX) 

PXM=INTGRL (MISSXO , MVELX) 

PYM=INTGRL (MISSYO , MVELY) 

C 

C SECOND ORDER PROP NAV MISSILE 

GDDOT = -20*GDOT + ( SOLOS-GAMMA) * 100 
GDOT = INTGRL (GDO, GDDOT) 

GAMMA = INTGRL (GO, GDOT) 

SOHDG = INTGRL (PHDG- 4*GDOT) 

SOXM = INTGRL (MISSXO, VM‘COS (SOHDG) ) 

SOYM = INTGRL(MISSYO,VM*SIN(SOHDG) ) 

C 

DYNAMIC 

C 

C THIRD ORDER PROP NAV MISSILE GEOMETRY UPDATE 

RM = ((XT-PXM)**2 + (YT-PYM) **2) “ . 5 

LOS = ATAN2 (YT-PYM,XT-PXM) 

RDOTM = TGTV*SIN (TGTHDG) /COS (LOS) 

VM*COS (PNHDG-LOS) 

C COMPUTE THE ERROR TERMS 

DELR(l) = RM - RKPl 
DELR(2) = RDOTM - RDKPl 
SDEL = LOS - SKPl 
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CALL KALMAN ( RNG , RM , RDOTM , RDDOTM , RK , RDK , RDDK , RKP 1 , 
+RDKP1 , RDDKPl , DELR , S , LOS , LOSD , LOSDD , SK , SDK , 

+SDDK , TIME , SKPl , SDKPl , SDDKPl , SDEL , GS , GR , RPHI , SPHI ) 

C 

BDDOT=SDDKP1+10* (SDK-BDOT) +33.33333* (LOS-B) 

U = -4*BDOT*RDOTM 



C 

C SECOND ORDER PROP NAV MISSILE GEOMETRY UPDATE 

SOR = ( (XT-SOXM) **2 + (YT-SOYM) * *2 ) * * . 5 
SOLOS = ATAN2 (YT-SOYM, XT-SOXM) 

SOU = -4*GDOT*SORD 
C 

IF (PXM .GT. XT) THEN 
CALL ENDRUN 
END IF 
C 

C STATEMENTS TO SAVE DATA FOR PLOTTING WITH DISSPLA 

IF (K .LE. 0) THEN 

WRITE (39,20) PXM , PYM , SOXM , SOYM 
WRITE (47,30) TIME, LOS, U 
WRITE (48,30) TIME , SOLOS , SOU 
K=10 

IF (RM .LT. .1*R) THEN 
K=3 
ENDIF 

NN=NN+1 



20 

30 

SAMPLE 


ENDIF 
K=K— 1 

FORMAT (4 (2X, FIO . 2) ) 

FORMAT (F5.2,2X,E11.3,2X,E14.6) 


C 


STATEMENTS TO SAVE DATA FOR GRAFAEL 


C 


SAVE 


(A) 0.1,XT,YT,PXM,PYM 


c 


SAVE 


(B) O.i, LOS,B,SK 


c 


SAVE 


(C 0.1, LOSD, BDOT, SDK 


c 


SAVE 


(D) 0.1, SDEL 


c 


SAVE 


E 0.1,RM,RK, RDDOTM, RDDK 


c 


SAVE 


(F) 0.1,U, BDDOT 


c 


SAVE 


G) 0 . 1 ,GS1 ,GS2 ,GS3 


c 

c 


SAVE 


(H) 0.1,GR11,GR12,GR21,GR22 



0 . 1 , RM , RKPl , RK , RDOTM , RDKPl , RDK , RDDOTM , RDDKPl , RDDK 
C 

0.1, PNHDG , LOS , B , SK , LOSD , BDOT , SDK , LOSDD , BDDOT , SDDK 
CONTROL FINTIM=10.0,DELT=.01 
C 

TERMINAL 



40 

C 

C 

(SC= 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



READ (1,40) MM 
WRITE (1,40) NK 
FORMAT (F6.1) 

STATEMENTS FOR PLOTTING USING GRAFAEL 

GRAPH (A/A,DE=TEK618) 

1600,LO=0.0) ,YT(SC=500,PO=16000) 

GRAPH (A/A,OV) PXM ( SC=1600 , AX=OMIT ) , PYM ( SC=500 ) 
B/B , DE=TEK618 ) TIME , LOS ( SC= . 025 , LO=- . 1 ) 
(B/B,OV) TIME(AX=OMIT) , B (PO=7 . 5 , SC= . 025 , LO=- 
(B/B,OV) TIME (AX=OMIT) , SK ( AX=OMIT, SC= . 025 , LO^ 
(C/C,DE=TEK618) TIME, LOSD 
IC/C.OV) TIME(AX=OMIT) , BDOT ( PO=7 . 5 ) 

(C/C,OV) TIME(AX=OMIT) , SDK ( AX=OMIT) 
(D/D,DE=TEK618) TIME, SDEL 

(E/E,DE=TEK618) TIME , RM ( SC=2000 . 0 , LO=0 . 0 ) 
E/E,OV) TIME (AX=OMIT) ,RK(SC=2000.0,LO=0.0) 

( F/E , DE=TEK618 ) TIME , RDDOTM ( SC=500 . 0 ) 

F/E,OV) TIME(AX=OMIT) , RDDK ( AX=OMIT) 
(G/F,DE=TEK618) TIME,U 

TIME,GS1 
TIME,GS2 
TIME,GS3 
TIME,GR11 



GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 

GRAPH 



(H/G,DE=TEK618 

(I/G,DE=TEK618 

J/G,DE=TEK618 

(K/H,DE=TEK618 



PRINT 

PRINT 



XT 



1 ) 

-. 1 ) 
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C GRAPH (L/H,DE=TEK618) TIME,GR12 

C GRAPH {M/H,DE=TEK618 TIME.GR21 

C GRAPH {N/H,DE=TEK618) TIME.GR22 

C GRAPH (0/F,DE=TEK618) TIME,BDDOT 

END 

STOP 

FORTRAN 

SUBROUTINE 

KALMAN { RNG , RM , RDOTM , RDDOTM , RK , RDK , RDDK , RKP 1 , RDKP 1 , 

* 

RDDKPl , DELR , S , LOS , LOSD , LOSDD , SK , SDK , SDDK , TIME , SKPl , SDKPl , 

* SDDKP1,SDEL,GS,GR,RPHI,SPHI) 

C SUBROUTINE TO ITERATE A KALMAN FILTER FOR RANGE 

VARIABLES 

C GIVEN THE COVARIANCE MATRIX, OBSERVATIONS 

REAL*8 

RNG { 3 ) , RM , RDOTM , RDDOTM , RK , RDK , RDDK , RKP 1 , RDKP 1 , RDDKP 1 , 

* 

DELR (2) ,S (3) , LOS, LOSD, RPHI (3,3) , A , B , C , D , TEMPI ( 3 ) ,SPHI (3,3) , 

* GR(3,3) ,GS(3) 

C 

C FIND THE NEW VALUES OF RPHI FROM THE PREVIOUS VALUES 

OF SIGMA 
C MATRIX 

C 

DO 110 N = 1.3 
TEMPI (N) = 0.0 
lie CONTINUE 

C CONSTANT GAIN INPUTS 

DO 121 L=l,3 
DO 120 M=l,2 

TEMPKL) = TEMPKL) + GR(L,M) * DELR(M) 

120 CONTINUE 

12 _ CONTINUE 

DO 125 N=l,3 

RNG(N) = TEMPI (N) + RNG(N) 

125 CONTINUE 

C SAVE THE VALUES OF RANGE MATRIX AT STEP K 

R:: = RNG(l) 

ROK = RNG 2) 

RODK = RNG (3) 

C 

C ZERO THE OLD RANGE TEMPORARY MATRIX 

C FIND THE ESTIMATE OF THE STEP K+1 FOR THE RANGE MATRIX 

DC 131 L = 1,3 
D^ 130 M = 1,3 

TEMPKL) = TEMPKL) + RPHI(L,M) * RNG(M) 

130 CONTINUE 

131 CONTINUE 
C 

C SAVE THE VALUES OF RNG(K+1/K) 

DO 132 N=l,3 

RNG (N) = TEMPI (N) 

132 CONTINUE 
C 

RKPl = RNG(l) 

RDKPl = RNG (2 

RDDKPl = RNG (3) 

C 

C SUBROUTINE TO ITERATE A KALMAN FILTER FOR SIGMA 

VARIABLES 

C GIVE!: THE COVARIANCE MATRIX, OBSERVATION. 

C WHERE H= (1 0 0) 

C 

C NOW FIND THE CURRENT VALUES OF THE SIGMA MATRIX 

S (1) = S (1) + GS (1) *SDEL 
S(2 = S(2 + GS2) *SDEL 
S\3) = S(3) + GS(3) *SDEL 
C 
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C STORE THE SIGMA MATRIX FOR USE IN THE PROGRAM 

SK = S(l) 

SDK = S(2) 

SDDK = S(3) 

C FIND THE NEXT VALUES OF THE SIGMA MATRIX 

C 

C S(K+1) = SPHI * S(K) 

C 

C ZERO A TEMPORARY MATRIX 

C 

DO 140 L=l,3 

TEMPI (L) = 0 

140 CONTINUE 

DO 142 L=l,3 

DO 141 M=l,3 

TEMPI (L)=TEMP1 (L) + SPHI(L,M) * S(M) 

141 CONTINUE 

142 CONTINUE 
C 

C INPUT BACK INTO SIGMA MATRIX 

DO 143 N = 1.3 
S (N) = TEMPI (N) 

143 CONTINUE 
C 

C STORE THE VALUE OF THE S MATRIX 

SKPl = S(l) 

SDKPl = S 2 

SDDKPl = S(3) 

C 

RETURN 

END 

SUBROUTINE ZERO (A, N) 

C CLEAR A TEMPORARY MATRIX 

C 

REAL* 8 A (3, 3) 

DO 201 L=1,N 

DO 200 M=1,N 
A (L,M) =0 

200 CONTINUE 

201 CONTINUE 
RETURN 

END 
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